DATA
Pop sive over time
dataset
data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
head(data)
## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5
## 1 no
## 2 no
## 3 no
## 4 no
## 5 no
## 6 no
tail(data)
## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Med 5 6 8/25/21
## 572 Med 5 6 8/25/21
## 573 Med 5 6 8/25/21
## 574 Med 5 6 8/25/21
## 575 Med 5 6 8/25/21
## 576 Med 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5
## 571 no no yes
## 572 extinct yes yes
## 573 no no no
## 574 extinct yes yes
## 575 no no no
## 576 no no no
dim(data)
## [1] 576 23
summary(data)
## Block Old_Label Label Population
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## Length:576 Min. :0.0 Min. :1.0 Length:576
## Class :character 1st Qu.:1.0 1st Qu.:2.0 Class :character
## Mode :character Median :2.5 Median :3.5 Mode :character
## Mean :2.5 Mean :3.5
## 3rd Qu.:4.0 3rd Qu.:5.0
## Max. :5.0 Max. :6.0
##
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## MC_Box1 MC_Box2 MC_Total_Adults HC_Box1
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 26.52 1st Qu.: 25.99 1st Qu.: 51.05 1st Qu.: 22.00
## Median : 58.87 Median : 54.29 Median :101.50 Median : 56.50
## Mean : 75.88 Mean : 70.57 Mean :128.37 Mean : 72.06
## 3rd Qu.:106.93 3rd Qu.: 99.40 3rd Qu.:170.25 3rd Qu.:101.00
## Max. :327.40 Max. :299.00 Max. :514.40 Max. :288.00
## NA's :142 NA's :141 NA's :5 NA's :144
## HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## Min. : 0.00 Min. : 0.0 Min. : 0.00 Min. :0.0000
## 1st Qu.: 25.00 1st Qu.: 62.0 1st Qu.: 67.25 1st Qu.:0.4183
## Median : 52.00 Median :100.0 Median :100.00 Median :0.9386
## Mean : 68.08 Mean :131.8 Mean :138.30 Mean :1.4482
## 3rd Qu.: 96.75 3rd Qu.:169.0 3rd Qu.:188.00 3rd Qu.:2.3587
## Max. :290.00 Max. :499.0 Max. :499.00 Max. :6.8571
## NA's :142 NA's :43 NA's :122 NA's :142
## First_Throw Extinction_finalstatus Less_than_5
## Length:576 Length:576 Length:576
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
#Remove populations killed due to the pathogen
data_status <- data.frame(data$Population,data$Extinction_finalstatus)
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)
## [1] 72 23
data<-data[data$Extinction_finalstatus!="kill",]
#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove
#Check variables
str(data)
## 'data.frame': 504 obs. of 26 variables:
## $ Block : chr "Block3" "Block3" "Block3" "Block3" ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : chr "High1" "High2" "High3" "High4" ...
## $ Genetic_Diversity : chr "High" "High" "High" "High" ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: chr "no" "no" "no" "no" ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)
#Order levels
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)
## [1] "High" "Medium" "Low"
data<-droplevels(data) #remove previous levels
str(data)
## 'data.frame': 504 obs. of 26 variables:
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
## $ Genetic_Diversity : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
head(data)
## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5 Nb_adults Lambda obs
## 1 no 100 NA 1
## 2 no 100 NA 2
## 3 no 100 NA 3
## 4 no 100 NA 4
## 5 no 100 NA 5
## 6 no 100 NA 6
tail(data)
## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Medium 5 6 8/25/21
## 572 Medium 5 6 8/25/21
## 573 Medium 5 6 8/25/21
## 574 Medium 5 6 8/25/21
## 575 Medium 5 6 8/25/21
## 576 Medium 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 571 no no yes 4 0.500000 499
## 572 extinct yes yes NA NA 500
## 573 no no no 46 1.314286 501
## 574 extinct yes yes NA NA 502
## 575 no no no 56 1.400000 503
## 576 no no no 133 3.243902 504
dim(data)
## [1] 504 26
#Add variable
data$gen_square <- data$Generation_Eggs^2
Heterozygosity remain
over time
#Upload He dataset
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]
#Merge dataset: add heterozygosity data to oridinal data
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)
### Creation of 4 He:
# 1- He_start: He was remaining for each population when we started the experiment
# 2- He_lost_gen_t: He was lost only during this specfic generation during the experiment (considering He=1 before this generation)
# 3- He_remain: He was remaining after each generation (He_remain[Gen=6]=He_end)
# 4- He_exp: He was lost during the entire experiment (considering He=1 at the beginning of the exp)
# 5- He_end: He was remaining for each population when we finished the experiment (considering the real He at the beginning of the exp, He_start, not 1)
#####
##### 1- He_start
#####
colnames(data)[28] <- "He_start"
#####
##### 2- He_lost_gen_t
#####
#Add He lost each generation
data$He_lost_gen_t <- 1 - (1/(2*data$Nb_adults))
#To consider extinct populations
is.na(data$He_lost_gen_t) <- sapply (data$He_lost_gen_t, is.infinite)
#####
##### 3- He_remain
#####
#Create new variable He remain per generation
data$He_remain <- NA
#Initiate with gen=1: He_start * He_lost_gen1
data$He_remain[data$Generation_Eggs=="1"] <- data$He_start[data$Generation_Eggs=="1"]*
data$He_lost_gen_t[data$Generation_Eggs=="1"]
#Compute for all the other genrations, except the first one
for(i in levels(data$Population)){
ifelse(sum(data$Generation_Eggs[data$Population==i])!=21,print("ERROR"),print(i))
for(j in 2:max(data$Generation_Eggs)){
data$He_remain[data$Population==i&data$Generation_Eggs==j] <-
data$He_remain[data$Population==i&data$Generation_Eggs==j-1]*
data$He_lost_gen_t[data$Population==i&data$Generation_Eggs==j]
}
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
#####
##### 4- He_exp
#####
# He at the end of the experiment
data$He_exp <- NA
#Compute the HE at the end of the experiment
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Eggs)!=21,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_adults))
#To consider extinct add this row:
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
#Compute for the whole experiment
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data$He_exp[data$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
#####
##### 5- He_end
#####
# Remaining He at the end of the experiment
data$He_end <- data$He_start * data$He_exp
#Save these values for Phenotyping
data_he <- merge(x = data_he,
y = data[data$Generation_Eggs==1,
c("Population","He_start",
"He_exp","He_end")],
by = "Population", all.x=TRUE)
Phenotyping
dataset
#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)
#Factors variables
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <- plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity,
levels = c("High", "Medium", "Low"))
#Add sum total
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx
data_phenotyping_4week <- data_phenotyping[data_phenotyping$Week=="4-week",]
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]
#Merge with He dataset
data_phenotyping <- merge(x = data_phenotyping,
y = data_he, by = "Population", all.x=TRUE)
data_phenotyping_4week <- merge(x = data_phenotyping_4week,
y = data_he, by = "Population", all.x=TRUE)
#Clean if less than 30 parents
pop_possible<-unique(data$Population[data$Generation_Eggs==5&data$Nb_adults>=30])
data_phenotyping <- data_phenotyping[data_phenotyping$Population %in% pop_possible,]
pop_possible<-unique(data$Population[data$Generation_Eggs==4&data$Nb_adults>=30])
data_phenotyping_4week <- data_phenotyping_4week[data_phenotyping_4week$Population %in% pop_possible,]
data_phenotyping_all <- rbind(data_phenotyping,data_phenotyping_4week)
Old He remain
# #Upload growth rate end of experiment
# data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
# data_he$Population <- as.factor(data_he$Population)
# data_he <- data_he[,c(1,3)]
#
#
# ### Additional: He remaining
# # New vector for the lost during exp
# data_he$He_remain_exp <- NA
#
# # He lost during exp
# for(i in levels(data$Population)){
# temp_data_pop <- subset(data, Population==i)
# ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
# temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))
#
# #To consider extinct add this row:
# is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
#
# temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
# data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
# rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
# }
#
# # Remaining He at the end of the experiment
# data_he$He_end <- data_he$He_remain * data_he$He_remain_exp
#
# # Remove kill population
# data_he <- data_he[!is.na(data_he$He_remain_exp),]
# data_he <- droplevels(data_he)
#
# ## Merge dataset He
# data_old_merged <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)
Survival data
library(survival)
library(lubridate)
library(ggsurvfit)
library("gtsummary")
library("tidycmprsk")
library("condSURV")
#########
data_survival <- data[data$Generation_Parents==0,c("Population","Block","Genetic_Diversity")]
data_survival$Gen_eggs_extinct <- max(data$Generation_Eggs)
data_survival$Survival <- "yes"
#Format data
for (i in unique(data_survival$Population)) {
if (data$Extinction_finalstatus[data$Population==i&data$Generation_Eggs==1] == "yes" ) {
data_survival$Gen_eggs_extinct[data_survival$Population==i]<-
min(data$Generation_Eggs[data$Population==i&data$First_Throw=="extinct"])
data_survival$Survival[data_survival$Population==i] <- "extinct"
}else{
print(i)
}
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low29"
## [1] "Low3"
## [1] "Low34"
## [1] "Low35"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med15"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med21"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
str(data_survival)
## 'data.frame': 84 obs. of 5 variables:
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 2 2 2 2 2 2 2 1 2 ...
## $ Genetic_Diversity: Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gen_eggs_extinct : int 6 6 6 6 6 6 6 6 6 6 ...
## $ Survival : chr "yes" "yes" "yes" "yes" ...
data_survival$Survival <- as.factor(data_survival$Survival)
#Note: the Surv() function in the {survival} package accepts by default TRUE/FALSE, where TRUE is event and FALSE is censored; 1/0 where 1 is event and 0 is censored; or 2/1 where 2 is event and 1 is censored. Please take care to ensure the event indicator is properly formatted.
data_survival$status <- ifelse(data_survival$Survival=="extinct", 1, 0)
PLOT
Population size
PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Population,
colour = Genetic_Diversity)) +
geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
SUM_Popsize <- Rmisc::summarySE(data,
measurevar = c("Nb_adults"),
groupvars = c("Generation_Eggs",
"Genetic_Diversity"),
na.rm=TRUE)
SUM_Popsize
## Generation_Eggs Genetic_Diversity N Nb_adults sd se ci
## 1 1 High 27 100.00000 0.00000 0.000000 0.00000
## 2 1 Medium 23 100.00000 0.00000 0.000000 0.00000
## 3 1 Low 34 100.00000 0.00000 0.000000 0.00000
## 4 2 High 27 344.29630 73.77294 14.197610 29.18360
## 5 2 Medium 23 282.04348 100.75735 21.009360 43.57075
## 6 2 Low 34 282.61765 78.53006 13.467794 27.40043
## 7 3 High 27 151.85185 80.96899 15.582489 32.03027
## 8 3 Medium 23 118.73913 88.54491 18.462891 38.28969
## 9 3 Low 34 104.61765 85.13341 14.600260 29.70445
## 10 4 High 26 114.00000 80.20224 15.728954 32.39439
## 11 4 Medium 23 62.65217 64.50483 13.450188 27.89398
## 12 4 Low 34 46.38235 39.63241 6.796903 13.82840
## 13 5 High 27 103.62963 51.56486 9.923661 20.39838
## 14 5 Medium 21 58.57143 51.11122 11.153383 23.26555
## 15 5 Low 31 51.51613 46.13377 8.285870 16.92200
## 16 6 High 27 147.66667 42.60282 8.198916 16.85311
## 17 6 Medium 19 69.73684 46.02396 10.558620 22.18284
## 18 6 Low 28 77.03571 46.14428 8.720450 17.89289
# SUM_Popsize <- Rmisc::summarySE(data[data$Extinction_finalstatus=="no",],
# measurevar = c("Nb_adults"),
# groupvars = c("Generation_Eggs",
# "Genetic_Diversity"),
# na.rm=TRUE)
# SUM_Popsize
PLOT_Pop_size_average <- PLOT_Pop_size +
geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 5, position = position_dodge(0.4)) +
geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
ymin = Nb_adults-ci, ymax = Nb_adults+ci,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 1, width=.2, position = position_dodge(0.4)) +
geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
linetype = "dashed", size = 1.5, position = position_dodge(0.4))
PLOT_Pop_size_average
# cowplot::save_plot(file = here::here("figures","PopulationSize_average.pdf"), PLOT_Pop_size_average, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)
#Prediction with models
# Test effect
# data$gen_square <- data$Generation_Eggs^2
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity +
gen_square*Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",],
family = "poisson")
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(2,6, length = 100),each=3),
Block = "Block4",
gen_square = (rep(seq(2,6, length = 100),each=3))^2,
Estimates = NA)
filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)
#Predcit estimate minimun values
tapply(filldata$Estimates, filldata$Genetic_Diversity, min)
## High Low Medium
## 85.85198 44.16308 51.99642
filldata[filldata$Estimates>=44.16&filldata$Estimates<=44.17,]
## Genetic_Diversity Generation_Eggs Block gen_square Estimates
## 198 Low 4.626263 Block4 21.40231 44.16308
filldata[filldata$Estimates>=51.99&filldata$Estimates<=52.00,]
## Genetic_Diversity Generation_Eggs Block gen_square Estimates
## 230 Medium 5.070707 Block4 25.71207 51.99642
filldata[filldata$Estimates>=85.85&filldata$Estimates<=85.86,]
## Genetic_Diversity Generation_Eggs Block gen_square Estimates
## 181 High 4.424242 Block4 19.57392 85.85198
## Change name vector
vector_names <- c(`Low` = "Strong bottleneck",
`Medium` = "Intermediate bottleneck",
`High` = "Diverse")
PLOT_Pop_size_predict<-ggplot2::ggplot(data[data$Extinction_finalstatus=="no",],
aes(x = factor(Generation_Eggs),
y = Nb_adults)) +
geom_point(aes(group = Population,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 0.7, alpha = 0.5) +
geom_line(aes(group = Population,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 0.05, alpha = 0.5) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.4) +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07", "#E7B800")) +
ylab("Population size") +
xlab("Generation") +
theme_LO
PLOT_Pop_size_predict
#
# cowplot::save_plot(file = here::here("figures","PopulationSize_predict.pdf"), PLOT_Pop_size_predict, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)
#
#
Extinction
#Percent of extinction
length(unique(data$Population[data$Genetic_Diversity=="Low"&
data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="Low"]))
## [1] 0.2352941
length(unique(data$Population[data$Genetic_Diversity=="Medium"&
data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="Medium"]))
## [1] 0.2173913
length(unique(data$Population[data$Genetic_Diversity=="High"&
data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="High"]))
## [1] 0
#Create data with percent
vector_names <- c(`Low` = "Strong bottleneck",
`Medium` = "Intermediate bottleneck",
`High` = "Diverse")
f_labels = data.frame(Genetic_Diversity = c("High","Medium","Low"),
label = c("Extinction = 0 %", "Extinction = 21.7 %", "Extinction = 23.5 %"))
f_labels$Genetic_Diversity <- factor(f_labels$Genetic_Diversity,
levels = c("High", "Medium", "Low"))
## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults)) +
facet_wrap(~ Genetic_Diversity, ncol=3, labeller = ggplot2::as_labeller(vector_names)) +
#geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
geom_line(aes(group = Population,
colour = Extinction_finalstatus),
position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
geom_text(x = 5, y = 420,
aes(label = label),
data = f_labels,
col="red",
size = 3,
vjust = 0) +
scale_color_manual(values = c("black", "red")) +
ylab("Population size") +
xlab("Generation") +
theme(legend.position = "none") +
theme_LO
PLOT_Extinction
#
#
# cowplot::save_plot(here::here("figures","Extinction.pdf"), PLOT_Extinction, base_height = 8/cm(1),
# base_width = 20/cm(1), dpi = 610)
# #
#
PLOT_Extinction
Growth rate G2
tapply(data[data$Generation_Eggs==2,]$Lambda,data[data$Generation_Eggs==2,]$Genetic_Diversity,mean)
## High Medium Low
## 3.442963 2.820435 2.826176
PLOT_Growth<-ggplot2::ggplot(data[data$Generation_Eggs==2,], aes(x = Genetic_Diversity, y = Lambda,
colour = Genetic_Diversity)) +
geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Genetic diversity") +
theme_LO
PLOT_Growth
#
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"), PLOT_Growth, base_height = 10/cm(1),
# base_width = 8/cm(1), dpi = 610)
Life stage
proportion
SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_adults"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_adults
## Genetic_Diversity Week N p_adults sd se ci
## 1 High 4-week 40 0.6315807 0.18084723 0.02859446 0.05783775
## 2 High 5-week 94 0.9646852 0.05234098 0.00539856 0.01072047
## 3 Medium 4-week 18 0.7261684 0.16721535 0.03941304 0.08315424
## 4 Medium 5-week 38 0.9594668 0.06046606 0.00980889 0.01987470
## 5 Low 4-week 27 0.6649974 0.19731065 0.03797245 0.07805350
## 6 Low 5-week 52 0.9538157 0.08241700 0.01142918 0.02294504
SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_pupae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_pupae
## Genetic_Diversity Week N p_pupae sd se ci
## 1 High 4-week 40 0.31153108 0.15526911 0.024550202 0.049657471
## 2 High 5-week 94 0.01417884 0.02356334 0.002430373 0.004826238
## 3 Medium 4-week 18 0.23410139 0.15149574 0.035707888 0.075337059
## 4 Medium 5-week 38 0.01931921 0.02904804 0.004712215 0.009547854
## 5 Low 4-week 27 0.28241476 0.17824104 0.034302504 0.070509807
## 6 Low 5-week 52 0.01796037 0.02546214 0.003530964 0.007088705
SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_larvae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_larvae
## Genetic_Diversity Week N p_larvae sd se ci
## 1 High 4-week 40 0.05688822 0.04383041 0.006930197 0.014017646
## 2 High 5-week 94 0.02113595 0.04275336 0.004409672 0.008756736
## 3 Medium 4-week 18 0.03973025 0.03784416 0.008919954 0.018819458
## 4 Medium 5-week 38 0.02121399 0.04478738 0.007265472 0.014721244
## 5 Low 4-week 27 0.05258789 0.08541162 0.016437473 0.033787710
## 6 Low 5-week 52 0.02822390 0.07073457 0.009809120 0.019692631
SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
p_pupae=SUM_Prop_pupae$p_pupae,
p_larvae=SUM_Prop_larvae$p_larvae)
SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop
## Genetic_Diversity Week N Stage Proportion
## 1 High 4-week 40 p_adults 0.63158069
## 2 High 5-week 94 p_adults 0.96468521
## 3 Medium 4-week 18 p_adults 0.72616835
## 4 Medium 5-week 38 p_adults 0.95946680
## 5 Low 4-week 27 p_adults 0.66499735
## 6 Low 5-week 52 p_adults 0.95381574
## 7 High 4-week 40 p_pupae 0.31153108
## 8 High 5-week 94 p_pupae 0.01417884
## 9 Medium 4-week 18 p_pupae 0.23410139
## 10 Medium 5-week 38 p_pupae 0.01931921
## 11 Low 4-week 27 p_pupae 0.28241476
## 12 Low 5-week 52 p_pupae 0.01796037
## 13 High 4-week 40 p_larvae 0.05688822
## 14 High 5-week 94 p_larvae 0.02113595
## 15 Medium 4-week 18 p_larvae 0.03973025
## 16 Medium 5-week 38 p_larvae 0.02121399
## 17 Low 4-week 27 p_larvae 0.05258789
## 18 Low 5-week 52 p_larvae 0.02822390
SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))
# New facet label names for supp variable
supp.weeks <- c("4 week", "5 week")
names(supp.weeks) <- c("4-week", "5-week")
PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity,
y = Proportion, fill = Stage)) +
facet_wrap(~ Week, ncol=2,
labeller = labeller(Week = supp.weeks)) +
geom_bar(stat="identity") +
xlab("Genetic history") +
labs(color="Stage of individuals") +
ylab("Proportion of each stage") +
scale_fill_manual(values= c("#619CFF","#F8766D","#00BA38"),
breaks = c("p_adults", "p_pupae", "p_larvae"),
labels = c( "Adult", "Pupae","Larvae")) +
ggplot2::scale_x_discrete(breaks=c("High", "Medium", "Low"),
labels=c("Diverse",
"Intermediate\nbottleneck",
"Strong \nbottleneck")) +
theme_LO +
theme(strip.background = element_rect(color="grey30",
fill="white", size=0.5, linetype="solid"),
strip.text.x = element_text(size = 12, color = "black", face = "bold"))
PLOT_prop
#
# cowplot::save_plot(here::here("figures","Life_stage_proportion.pdf"), PLOT_prop,
# base_height = 10/cm(1), base_width = 18/cm(1), dpi = 610)
Growth rate and
He
# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_phenotyping,
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity",
"He_end",
"Population"),
na.rm=TRUE)
PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_end,
y = Lambda,
colour = Genetic_Diversity,
size = N)) +
geom_point(alpha = 0.8) +
ylab("Growth rate") +
xlab("Expected heterozygosity") +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07","#E7B800")) +
labs(size="Replicates") +
theme_LO
PLOT_He_Final
#
#
# cowplot::save_plot(here::here("figures","Heterozygosity_Growthrate.pdf"),
# PLOT_He_Final,
# base_height = 12/cm(1), base_width = 16/cm(1), dpi = 610)
#
#
# # ALL WITHOUT PSEUDO
# SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
# measurevar = c("Lambda"),
# groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_end",
# "Population"),
# na.rm=TRUE)
#
# SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_end)
#
# SUM_POP_He_ALL <- SUM_POP_He_ALL[!is.na(SUM_POP_He_ALL$He_end),]
# SUM_POP_He_ALL <- SUM_POP_He_ALL[SUM_POP_He_ALL$HE!="-Inf",]
#
#
# facet_names <- c('Initial'="Before rescue experiment",'Final'="After rescue experiment")
#
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
# colour = Genetic_Diversity,
# shape = Phenotyping)) +
# #facet_wrap(~ Phenotyping, ncol=1) +
# geom_point(size = 3, alpha = 0.8) +
# scale_shape_manual(values = c(1, 19)) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL
#
#
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
# colour = Genetic_Diversity,
# shape = Phenotyping)) +
# facet_wrap(~ Phenotyping, ncol=1) +
# geom_point(size = 3, alpha = 0.8) +
# scale_shape_manual(values = c(1, 19)) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL
#
#
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
# colour = Genetic_Diversity)) +
# facet_wrap(~ Phenotyping, ncol=1) +
# geom_point(size = 3, alpha = 0.8) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL
Lifestage and He
SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_adults"),
groupvars = c("Genetic_Diversity","Week", "Population", "He_end"),
na.rm = TRUE)
colnames(SUM_Prop_adults)[6]<-"prop"
SUM_Prop_adults$lifestage <- "adults"
SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_pupae"),
groupvars = c("Genetic_Diversity","Week", "Population", "He_end"),
na.rm = TRUE)
colnames(SUM_Prop_pupae)[6]<-"prop"
SUM_Prop_pupae$lifestage <- "pupae"
SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_larvae"),
groupvars = c("Genetic_Diversity","Week",
"Population", "He_end"),
na.rm = TRUE)
colnames(SUM_Prop_larvae)[6]<-"prop"
SUM_Prop_larvae$lifestage <- "larvae"
SUM_Prop_POP <- rbind(SUM_Prop_adults, SUM_Prop_pupae, SUM_Prop_larvae)
PLOT_Proportion_lifestage <- ggplot2::ggplot(SUM_Prop_POP, aes(x = He_end,
y = prop,
colour = Genetic_Diversity,
size = N)) +
facet_grid( Week ~ factor(lifestage, levels=c('larvae','pupae','adults'))) +
geom_point(alpha = 0.8) +
ylab("Proportion") +
xlab("Expected heterozygosity") +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07", "#E7B800")) +
labs(size="Replicates") +
theme_LO
PLOT_Proportion_lifestage
#
# cowplot::save_plot(here::here("figures","Life_stage_proportion_He.pdf"),
# PLOT_Proportion_lifestage,
# base_height = 14/cm(1), base_width = 36/cm(1), dpi = 610)
#
#
He over time
#Compute for all the other genrations, except the first one
data$ID_extinction <- "extant"
for(i in unique(data$Population[data$Extinction_finalstatus=="yes"])){
gen_all <- data$Generation_Eggs[data$Population==i&!is.na(data$He_remain)]
maxgen <- max(gen_all)
data$ID_extinction[data$Population==i&data$Generation_Eggs==maxgen] <- "willextinct"
}
data[data$ID_extinction=="willextinct",]
## Population Block Old_Label Label
## 205 Low16 Block4 B4-O1 7/14 BE B4 | G5 Low 16
## 273 Low27 Block5 B5-B1 6/16 BE B5 | G4 Low 27
## 278 Low28 Block5 B5_E1 5/12 BE B5 | G3 Low 28
## 299 Low30 Block5 B5-D1 6/16 BE B5 | G4 Low 30
## 303 Low31 Block5 B(2)5-P1 6/16 BE B5 | G4 Low 31
## 310 Low32 Block5 B(2)5_Q1 5/12 BE B5 | G3 Low 32
## 317 Low33 Block5 B(2)5-T1 7/21 BE B5 | G5 Low 33
## 336 Low36 Block5 B(2)5-X1 6/16 BE B5 | G4 Low 36
## 402 Med14 Block4 B4-Backup Fam L 6/9 BE B4 | G4 Med 14
## 409 Med17 Block5 B5_Backup_Fam_F 5/12 BE B5 | G3 Med 17
## 437 Med20 Block5 B5 - Backup Fam I 6/16 BE B5 | G4 Med 20
## 449 Med22 Block5 B5_Backup_Fam_B 5/12 BE B5 | G3 Med 22
## 479 Med5 Block3 B3 - Backup Fam E 7/7 BE B3 | G5 Med 5
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 205 Low 4 5 7/14/21
## 273 Low 3 4 6/16/21
## 278 Low 2 3 5/12/21
## 299 Low 3 4 6/16/21
## 303 Low 3 4 6/16/21
## 310 Low 2 3 5/12/21
## 317 Low 4 5 7/21/21
## 336 Low 3 4 6/16/21
## 402 Medium 3 4 6/9/21
## 409 Medium 2 3 5/12/21
## 437 Medium 3 4 6/16/21
## 449 Medium 2 3 5/12/21
## 479 Medium 4 5 7/7/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 205 7/14/21 7/15/21 <NA> DSC_0994 0.000000 2.000000
## 273 6/16/21 6/17/21 DSC_0526 DSC_0527 5.090878 0.000000
## 278 5/12/21 5/13/21 DSC_0918 DSC_0919 3.085768 3.638847
## 299 6/16/21 6/17/21 DSC_0559 <NA> 1.383920 0.000000
## 303 6/16/21 6/17/21 <NA> DSC_0541 0.000000 1.000000
## 310 5/12/21 5/13/21 DSC_0953 DSC_0954 88.254170 70.952124
## 317 7/21/21 7/22/21 DSC_0117 DSC_0118 1.000000 1.000000
## 336 6/16/21 6/17/21 DSC_0540 <NA> 1.000000 0.000000
## 402 6/9/21 6/10/21 DSC_0376 DSC_0377 10.974990 1.000000
## 409 5/12/21 5/13/21 DSC_0916 DSC_0917 1.543122 1.084327
## 437 6/16/21 6/17/21 DSC_0542 <NA> 1.000000 0.000000
## 449 5/12/21 5/13/21 DSC_0922 DSC_0923 9.168447 6.251069
## 479 7/7/21 7/8/21 <NA> DSC_0830 NA NA
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 205 2.0 0 2 2 5 0.400000000
## 273 5.1 3 0 3 19 0.157894737
## 278 6.7 1 1 2 374 0.005347594
## 299 1.4 1 0 1 146 0.006849315
## 303 1.0 1 1 2 272 0.007352941
## 310 159.2 71 56 127 276 0.460144928
## 317 2.0 1 1 2 3 0.666666667
## 336 1.0 1 0 1 22 0.045454545
## 402 12.0 7 1 8 138 0.057971014
## 409 2.6 3 2 5 416 0.012019231
## 437 1.0 1 0 1 20 0.050000000
## 449 15.4 10 6 16 200 0.080000000
## 479 0.0 NA NA 1 11 0.090909091
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 205 no yes yes 2 0.400000000 378
## 273 no yes yes 3 0.157894737 319
## 278 no yes yes 2 0.005347594 236
## 299 no yes yes 1 0.006849315 322
## 303 no yes yes 2 0.007352941 323
## 310 no yes yes 127 0.460144928 240
## 317 no yes yes 2 0.666666667 409
## 336 no yes yes 1 0.045454545 328
## 402 no yes yes 8 0.057971014 308
## 409 no yes yes 5 0.012019231 245
## 437 no yes yes 1 0.050000000 332
## 449 no yes yes 16 0.080000000 250
## 479 no yes yes 1 0.090909091 359
## gen_square He_start He_lost_gen_t He_remain He_exp He_end
## 205 25 0.5544299 0.7500000 0.3575672 0.6449276 0.3575672
## 273 16 0.5367998 0.8333333 0.4324691 0.8056432 0.4324691
## 278 9 0.5386585 0.7500000 0.4014365 0.7452523 0.4014365
## 299 16 0.5540444 0.5000000 0.2742819 0.4950540 0.2742819
## 303 16 0.5532775 0.7500000 0.4116634 0.7440450 0.4116634
## 310 9 0.5528880 0.9960630 0.5469651 0.9892872 0.5469651
## 317 25 0.5523333 0.7500000 0.3359893 0.6083089 0.3359893
## 336 16 0.5427371 0.5000000 0.2635269 0.4855518 0.2635269
## 402 16 0.6802331 0.9375000 0.6315974 0.9285014 0.6315974
## 409 9 0.6825509 0.9000000 0.6104897 0.8944237 0.6104897
## 437 16 0.6819650 0.5000000 0.3302071 0.4841994 0.3302071
## 449 9 0.6809168 0.9687500 0.6546991 0.9614965 0.6546991
## 479 25 0.6807065 0.5000000 0.3209456 0.4714889 0.3209456
## ID_extinction
## 205 willextinct
## 273 willextinct
## 278 willextinct
## 299 willextinct
## 303 willextinct
## 310 willextinct
## 317 willextinct
## 336 willextinct
## 402 willextinct
## 409 willextinct
## 437 willextinct
## 449 willextinct
## 479 willextinct
dim(data[data$ID_extinction=="willextinct",])
## [1] 13 33
data[data$ID_extinction=="willextinct","Population"]
## [1] Low16 Low27 Low28 Low30 Low31 Low32 Low33 Low36 Med14 Med17 Med20 Med22
## [13] Med5
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9
data$Population[data$Genetic_Diversity=="Medium"&data$Extinction_finalstatus=="yes"]
## [1] Med14 Med14 Med14 Med14 Med14 Med14 Med17 Med17 Med17 Med17 Med17 Med17
## [13] Med20 Med20 Med20 Med20 Med20 Med20 Med22 Med22 Med22 Med22 Med22 Med22
## [25] Med5 Med5 Med5 Med5 Med5 Med5
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9
PLOT_He_extinction <-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = He_remain,
group = Population,
colour = Extinction_finalstatus,
shape = ID_extinction)) +
geom_line(position = position_dodge(0.5),
size = 0.25, alpha = 0.6) +
geom_point(position = position_dodge(0.5), size = 3, stroke = 0.6,
fill= "white", alpha = 0.8) +
ylab("Expected heterozygozity") +
scale_color_manual(values = c("black", "red")) +
scale_shape_manual(values = c(21, 13), guide=FALSE) +
labs(color="Extinct populations") +
xlab("Generation") +
theme_LO
PLOT_He_extinction
# #
# cowplot::save_plot(here::here("figures","Heterozygosity_over_time.pdf"),
# PLOT_He_extinction,
# base_height = 12/cm(1), base_width = 14/cm(1), dpi = 610)
#
Survival vs He
s2 <- survfit(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
summary(s2)
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity,
## data = data_survival)
##
## Genetic_Diversity=High
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
##
## Genetic_Diversity=Medium
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 23 2 0.913 0.0588 0.805 1.000
## 5 21 2 0.826 0.0790 0.685 0.996
## 6 19 1 0.783 0.0860 0.631 0.971
##
## Genetic_Diversity=Low
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 34 2 0.941 0.0404 0.865 1.000
## 5 32 4 0.824 0.0654 0.705 0.962
## 6 28 2 0.765 0.0727 0.635 0.921
# Create a summary dataset
SUM_survprob <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity), each =6),
Generation_Eggs = rep(seq(1:6), n=3),
SurvProb = 1,
se = 0)
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Medium"&
SUM_survprob$Generation_Eggs==4]<-summary(s2)$surv[1]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Medium"&
SUM_survprob$Generation_Eggs==5]<-summary(s2)$surv[2]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Medium"&
SUM_survprob$Generation_Eggs==6]<-summary(s2)$surv[3]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Low"&
SUM_survprob$Generation_Eggs==4]<-summary(s2)$surv[4]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Low"&
SUM_survprob$Generation_Eggs==5]<-summary(s2)$surv[5]
SUM_survprob$SurvProb[SUM_survprob$Genetic_Diversity=="Low"&
SUM_survprob$Generation_Eggs==6]<-summary(s2)$surv[6]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Medium"&
SUM_survprob$Generation_Eggs==4]<-summary(s2)$std.err[1]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Medium"&
SUM_survprob$Generation_Eggs==5]<-summary(s2)$std.err[2]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Medium"&
SUM_survprob$Generation_Eggs==6]<-summary(s2)$std.err[3]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Low"&
SUM_survprob$Generation_Eggs==4]<-summary(s2)$std.err[4]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Low"&
SUM_survprob$Generation_Eggs==5]<-summary(s2)$std.err[5]
SUM_survprob$se[SUM_survprob$Genetic_Diversity=="Low"&
SUM_survprob$Generation_Eggs==6]<-summary(s2)$std.err[6]
data_survprob_he <- merge(x = data, y = SUM_survprob, by = c("Genetic_Diversity", "Generation_Eggs"), all.x=TRUE)
head(data_survprob_he)
## Genetic_Diversity Generation_Eggs Population Block Old_Label
## 1 High 1 High1 Block3 B3_All_Red_Mix
## 2 High 1 High26 Block5 B5_All_Red_Mix
## 3 High 1 High14 Block4 B4_All_Red_Mix
## 4 High 1 High33 Block5 B5_All_Red_Mix
## 5 High 1 High4 Block3 B3_All_Red_Mix
## 6 High 1 High19 Block4 B4_All_Red_Mix
## Label Generation_Parents Date_Census Date_Start_Eggs
## 1 BE B3 2/17 | G1 High1 0 2/17/21 2/17/21
## 2 BE B5 3/3 | G1 High26 0 3/3/21 3/3/21
## 3 BE B4 2/24 | G1 High14 0 2/24/21 2/24/21
## 4 BE B5 3/3 | G1 High33 0 3/3/21 3/3/21
## 5 BE B3 2/17 | G1 High4 0 2/17/21 2/17/21
## 6 BE B4 2/24 | G1 High19 0 2/24/21 2/24/21
## Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1
## 1 2/18/21 DSC_0155 <NA> NA NA 101.9 NA
## 2 3/4/21 DSC_0675 <NA> NA NA 101.5 NA
## 3 2/25/21 DSC_0416 <NA> NA NA 107.5 NA
## 4 3/4/21 DSC_0682 <NA> NA NA 82.5 NA
## 5 2/18/21 DSC_0158 <NA> NA NA 101.9 NA
## 6 2/25/21 DSC_0422 <NA> NA NA 101.8 NA
## HC_Box2 HC_Total_Adults Nb_parents Growth_rate First_Throw
## 1 NA 100 NA NA no
## 2 NA 100 NA NA no
## 3 NA 100 NA NA no
## 4 NA 100 NA NA no
## 5 NA 100 NA NA no
## 6 NA 100 NA NA no
## Extinction_finalstatus Less_than_5 Nb_adults Lambda obs gen_square He_start
## 1 no no 100 NA 1 1 0.9986008
## 2 no no 100 NA 59 1 0.9986008
## 3 no no 100 NA 28 1 0.9986008
## 4 no no 100 NA 63 1 0.9986008
## 5 no no 100 NA 4 1 0.9986008
## 6 no no 100 NA 33 1 0.9986008
## He_lost_gen_t He_remain He_exp He_end ID_extinction SurvProb se
## 1 0.995 0.9936078 0.9679591 0.9666047 extant 1 0
## 2 0.995 0.9936078 0.9220546 0.9207645 extant 1 0
## 3 0.995 0.9936078 0.9766045 0.9752381 extant 1 0
## 4 0.995 0.9936078 0.9358965 0.9345870 extant 1 0
## 5 0.995 0.9936078 0.9812175 0.9798447 extant 1 0
## 6 0.995 0.9936078 0.9804728 0.9791010 extant 1 0
Plot_survival <- ggplot2::ggplot(data_survprob_he, aes(x = He_remain, y = SurvProb , colour = Genetic_Diversity)) +
geom_point(alpha = 0.8) +
ylab("Survival probability") +
xlab("Expected heterozygosity") +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07", "#E7B800")) +
labs(size="Replicates") +
theme_LO
#
#
# cowplot::save_plot(here::here("figures","Survival_function_He.pdf"),
# Plot_survival,
# base_height = 12/cm(1), base_width = 14/cm(1), dpi = 610)
#
ANALYSIS
Probability of
extinction
############
###### Clean dataset
############
#Prepare clean dataset
data_proba_extinction <- data[data$Generation_Eggs==6,
c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)
############
###### Analysis
############
#Analysis
mod0 <- glm(y ~ Genetic_Diversity + Block,
data = data_proba_extinction, family = binomial)
summary(mod0)
##
## Call:
## glm(formula = y ~ Genetic_Diversity + Block, family = binomial,
## data = data_proba_extinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26011 -0.44258 -0.30957 -0.00005 2.47474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.3144 1855.7435 -0.011 0.99084
## Genetic_DiversityMedium 18.3002 1855.7433 0.010 0.99213
## Genetic_DiversityLow 18.5062 1855.7432 0.010 0.99204
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 72.388 on 83 degrees of freedom
## Residual deviance: 46.833 on 79 degrees of freedom
## AIC: 56.833
##
## Number of Fisher Scoring iterations: 18
#Test convergence
#Test genetic diversity effect
mod1 <- glm(y ~ Block ,
data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 81 58.903
## 2 79 46.833 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#We keep genetic diversity
# Perform the lrtest
(A <- logLik(mod1))
## 'log Lik.' -29.45146 (df=3)
(B <- logLik(mod0))
## 'log Lik.' -23.41669 (df=5)
#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))
## [1] 12.06954
#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))
## [1] 0.002394043
lmtest::lrtest(mod1, mod0)
## Likelihood ratio test
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -29.451
## 2 5 -23.417 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,
data = data_proba_extinction, family = binomial)
summary(mod0)
##
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial,
## data = data_proba_extinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26011 -0.44258 -0.30957 -0.00005 2.47474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Genetic_DiversityHigh -21.3144 1855.7435 -0.011 0.99084
## Genetic_DiversityMedium -3.0142 1.1293 -2.669 0.00760 **
## Genetic_DiversityLow -2.8082 1.0659 -2.635 0.00843 **
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.449 on 84 degrees of freedom
## Residual deviance: 46.833 on 79 degrees of freedom
## AIC: 56.833
##
## Number of Fisher Scoring iterations: 18
############
###### Predicted value
############
#Extract probability of extinction
(p_high <- plogis(-21.3144))
## [1] 5.536989e-10
(p_med <- plogis(-3.0142))
## [1] 0.04678847
(p_low <- plogis(-2.8082))
## [1] 0.05688267
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
Block = "Block4")
data_predict_extinct$predict <- predict(mod0, type = "response",
re.form = NA,
newdata = data_predict_extinct)
data_predict_extinct
## Genetic_Diversity Block predict
## 1 High Block4 1.160723e-09
## 2 Medium Block4 9.329490e-02
## 3 Low Block4 1.122446e-01
p_high
## [1] 5.536989e-10
p_med
## [1] 0.04678847
p_low
## [1] 0.05688267
############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High -20.07 1855.743 Inf -4451.10 4410.961
## Medium -1.77 0.650 Inf -3.32 -0.216
## Low -1.56 0.532 Inf -2.83 -0.290
##
## Results are averaged over the levels of: Block
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium -18.300 1855.743 Inf -0.010 0.9999
## High - Low -18.506 1855.743 Inf -0.010 0.9999
## Medium - Low -0.206 0.751 Inf -0.274 0.9594
##
## Results are averaged over the levels of: Block
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
Time to
extinction
data_timeextinction<-data[data$First_Throw=="extinct",]
data_timeextinction <- data_timeextinction[complete.cases(data_timeextinction$Nb_adults), ]
data_timeextinction
## Population Block Old_Label Label
## 207 Low16 Block4 B4-O1 8/18 BE B4 | G6 Low 16
## 274 Low27 Block5 B5-B1 7/21 BE B5 | G5 Low 27
## 277 Low28 Block5 B5-E1 6/16 BE B5 | G4 Low 28
## 296 Low30 Block5 B5-D1 7/21 BE B5 | G5 Low 30
## 301 Low31 Block5 B(2)5-P1 7/21 BE B5 | G5 Low 31
## 309 Low32 Block5 B(2)5-Q1 6/16 BE B5 | G4 Low 32
## 318 Low33 Block5 B(2)5-T1 8/25 BE B5 | G6 Low 33
## 335 Low36 Block5 B(2)5-X1 7/21 BE B5 | G5 Low 36
## 397 Med14 Block4 B4-Backup Fam L 7/14 BE B4 | G5 Med 14
## 410 Med17 Block5 B5 - Backup Fam F 6/16 BE B5 | G4 Med 17
## 438 Med20 Block5 B5 - Backup Fam I 7/21 BE B5 | G5 Med 20
## 448 Med22 Block5 B5-Backup Fam B 6/16 BE B5 | G4 Med 22
## 476 Med5 Block3 B3 - Backup Fam E 8/11 BE B3 | G6 Med 5
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 207 Low 5 6 8/18/21
## 274 Low 4 5 7/21/21
## 277 Low 3 4 6/16/21
## 296 Low 4 5 7/21/21
## 301 Low 4 5 7/21/21
## 309 Low 3 4 6/16/21
## 318 Low 5 6 8/25/21
## 335 Low 4 5 7/21/21
## 397 Medium 4 5 7/14/21
## 410 Medium 3 4 6/16/21
## 438 Medium 4 5 7/21/21
## 448 Medium 3 4 6/16/21
## 476 Medium 5 6 8/11/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 207 8/18/21 8/19/21 <NA> <NA> NA NA
## 274 7/21/21 7/22/21 <NA> <NA> 0 0
## 277 6/16/21 6/17/21 <NA> <NA> 0 0
## 296 7/21/21 7/22/21 <NA> <NA> 0 0
## 301 7/21/21 7/22/21 <NA> <NA> 0 0
## 309 6/16/21 6/17/21 <NA> <NA> 0 0
## 318 8/25/21 8/26/21 <NA> <NA> NA NA
## 335 7/21/21 7/22/21 <NA> <NA> 0 0
## 397 7/14/21 7/15/21 <NA> <NA> 0 0
## 410 6/16/21 6/17/21 <NA> <NA> 0 0
## 438 7/21/21 7/22/21 <NA> <NA> 0 0
## 448 6/16/21 6/17/21 <NA> <NA> 0 0
## 476 8/11/21 8/12/21 <NA> <NA> 0 0
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 207 0 0 0 0 2 0
## 274 0 0 0 0 3 0
## 277 0 0 0 0 2 0
## 296 0 0 0 0 1 0
## 301 0 0 0 0 2 0
## 309 0 0 0 0 127 0
## 318 0 NA NA 0 2 0
## 335 0 0 0 0 1 0
## 397 0 NA NA 0 8 NA
## 410 0 0 0 0 5 0
## 438 0 0 0 0 1 0
## 448 0 0 0 0 16 0
## 476 0 NA NA 0 1 0
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 207 extinct yes yes 0 0 462
## 274 extinct yes yes 0 0 403
## 277 extinct yes yes 0 0 320
## 296 extinct yes yes 0 0 406
## 301 extinct yes yes 0 0 407
## 309 extinct yes yes 0 0 324
## 318 extinct yes yes 0 0 493
## 335 extinct yes yes 0 0 412
## 397 extinct yes yes 0 0 392
## 410 extinct yes yes 0 0 329
## 438 extinct yes yes 0 0 416
## 448 extinct yes yes 0 0 334
## 476 extinct yes yes 0 0 443
## gen_square He_start He_lost_gen_t He_remain He_exp He_end
## 207 36 0.5544299 NA NA 0.6449276 0.3575672
## 274 25 0.5367998 NA NA 0.8056432 0.4324691
## 277 16 0.5386585 NA NA 0.7452523 0.4014365
## 296 25 0.5540444 NA NA 0.4950540 0.2742819
## 301 25 0.5532775 NA NA 0.7440450 0.4116634
## 309 16 0.5528880 NA NA 0.9892872 0.5469651
## 318 36 0.5523333 NA NA 0.6083089 0.3359893
## 335 25 0.5427371 NA NA 0.4855518 0.2635269
## 397 25 0.6802331 NA NA 0.9285014 0.6315974
## 410 16 0.6825509 NA NA 0.8944237 0.6104897
## 438 25 0.6819650 NA NA 0.4841994 0.3302071
## 448 16 0.6809168 NA NA 0.9614965 0.6546991
## 476 36 0.6807065 NA NA 0.4714889 0.3209456
## ID_extinction
## 207 extant
## 274 extant
## 277 extant
## 296 extant
## 301 extant
## 309 extant
## 318 extant
## 335 extant
## 397 extant
## 410 extant
## 438 extant
## 448 extant
## 476 extant
tapply(data_timeextinction$Generation_Eggs,data_timeextinction$Genetic_Diversity,mean)
## High Medium Low
## NA 4.8 5.0
#Test time of extinction
mod1 <- glm(Generation_Eggs ~ Genetic_Diversity + Block,
data = data_timeextinction, family = "poisson")
mod0 <- glm(Generation_Eggs ~ Block,
data = data_timeextinction, family = "poisson")
anova(mod0, mod1, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ Genetic_Diversity + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 10 0.95459
## 2 9 0.74888 1 0.20571 0.6502
#We keep genetic diversity
# Perform the lrtest
(A <- logLik(mod1))
## 'log Lik.' -22.83384 (df=4)
(B <- logLik(mod0))
## 'log Lik.' -22.9367 (df=3)
#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))
## [1] -0.2057062
#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))
## [1] 1
lmtest::lrtest(mod0,mod1)
## Likelihood ratio test
##
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ Genetic_Diversity + Block
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -22.937
## 2 4 -22.834 1 0.2057 0.6502
Survival
analysis
str(data_survival)
## 'data.frame': 84 obs. of 6 variables:
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 2 2 2 2 2 2 2 1 2 ...
## $ Genetic_Diversity: Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gen_eggs_extinct : int 6 6 6 6 6 6 6 6 6 6 ...
## $ Survival : Factor w/ 2 levels "extinct","yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ status : num 0 0 0 0 0 0 0 0 0 0 ...
#Note: the Surv() function in the {survival} package accepts by default TRUE/FALSE, where TRUE is event and FALSE is censored; 1/0 where 1 is event and 0 is censored; or 2/1 where 2 is event and 1 is censored. Please take care to ensure the event indicator is properly formatted.
# The Surv() function from the {survival} package creates a survival object for use as the response in a model formula. There will be one entry for each subject that is the survival time, which is followed by a + if the subject was censored. Let’s look at the first 10 observations:
Surv(data_survival$Gen_eggs_extinct, data_survival$status)
## [1] 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+
## [26] 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 5 4 6+ 6+ 5
## [51] 5 4 6 6+ 6+ 5 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 6+ 5 6+ 4 6+ 6+ 6+ 5 6+ 4
## [76] 6+ 6+ 6+ 6+ 6 6+ 6+ 6+ 6+
#Surv(data_survival_all$Gen_eggs_extinct, data_survival$status)
# The survfit() function creates survival curves using the Kaplan-Meier method based on a formula. Let’s generate the overall survival curve for the entire cohort, assign it to object s1, and look at the structure using str():
s1 <- survfit(Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival)
str(s1)
## List of 16
## $ n : int 84
## $ time : num [1:3] 4 5 6
## $ n.risk : num [1:3] 84 80 74
## $ n.event : num [1:3] 4 6 3
## $ n.censor : num [1:3] 0 0 71
## $ surv : num [1:3] 0.952 0.881 0.845
## $ std.err : num [1:3] 0.0244 0.0401 0.0467
## $ cumhaz : num [1:3] 0.0476 0.1226 0.1632
## $ std.chaz : num [1:3] 0.0238 0.0388 0.0453
## $ type : chr "right"
## $ logse : logi TRUE
## $ conf.int : num 0.95
## $ conf.type: chr "log"
## $ lower : num [1:3] 0.908 0.814 0.771
## $ upper : num [1:3] 0.999 0.953 0.926
## $ call : language survfit(formula = Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival)
## - attr(*, "class")= chr "survfit"
# time: the timepoints at which the curve has a step, i.e. at least one event occurred
# surv: the estimate of survival at the corresponding time
# The {ggsurvfit} package works best if you create the survfit object using the included ggsurvfit::survfit2() function, which uses the same syntax to what we saw previously with survival::survfit(). The ggsurvfit::survfit2() tracks the environment from the function call, which allows the plot to have better default values for labeling and p-value reporting.
#
####
#### PLOT
####
ggsurvfit::survfit2(Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival) %>%
ggsurvfit() +
ggplot2::labs( x = "Generation",
y = "Overall survival probability") +
add_confidence_interval() +
add_risktable()
summary(survfit(Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival))
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ 1, data = data_survival)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 84 4 0.952 0.0232 0.908 0.999
## 5 80 6 0.881 0.0353 0.814 0.953
## 6 74 3 0.845 0.0395 0.771 0.926
##### COMPARISONS ACROSS GROUPS
# We can conduct between-group significance tests using a log-rank test. The log-rank test equally weights observations over the entire follow-up time and is the most common way to compare survival times between groups. There are versions that more heavily weight the early or late follow-up that could be more appropriate depending on the research question (see ?survdiff for different test options).
#
# We get the log-rank p-value using the survdiff function. For example, we can test whether there was a difference in survival time according to the demography history in the lung data:
#Comparison
survdiff(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
## Call:
## survdiff(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity,
## data = data_survival)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Genetic_Diversity=High 27 0 4.41 4.405 7.00
## Genetic_Diversity=Medium 23 5 3.44 0.707 1.01
## Genetic_Diversity=Low 34 8 5.15 1.571 2.73
##
## Chisq= 7 on 2 degrees of freedom, p= 0.03
summary(survfit(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival))
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity,
## data = data_survival)
##
## Genetic_Diversity=High
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
##
## Genetic_Diversity=Medium
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 23 2 0.913 0.0588 0.805 1.000
## 5 21 2 0.826 0.0790 0.685 0.996
## 6 19 1 0.783 0.0860 0.631 0.971
##
## Genetic_Diversity=Low
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 34 2 0.941 0.0404 0.865 1.000
## 5 32 4 0.824 0.0654 0.705 0.962
## 6 28 2 0.765 0.0727 0.635 0.921
####
#### PLOT
####
plot_survival <- ggsurvfit::survfit2(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival) %>%
ggsurvfit(size = 2) +
labs( x = "Generation",
y = "Overall survival probability") +
add_confidence_interval(alpha = 0.1) +
add_censor_mark(size = 8, shape = "x") +
scale_fill_manual(name="Demographic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07", "#E7B800")) +
scale_color_manual(name="Demographic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07", "#E7B800")) #+ add_risktable()
plot_survival
# cowplot::save_plot(file = here::here("figures","Extinction_Survival_analysis.pdf"),
# plot_survival, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)
# COX representation
coxph(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
## Call:
## coxph(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity,
## data = data_survival)
##
## coef exp(coef) se(coef) z p
## Genetic_DiversityMedium 1.967e+01 3.475e+08 7.239e+03 0.003 0.998
## Genetic_DiversityLow 1.974e+01 3.727e+08 7.239e+03 0.003 0.998
##
## Likelihood ratio test=11.1 on 2 df, p=0.003883
## n= 84, number of events= 13
coxph(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival) %>%
gtsummary::tbl_regression(exp = TRUE)
| Characteristic |
HR |
95% CI |
p-value |
| Genetic_Diversity |
|
|
|
| High |
— |
— |
|
| Medium |
347,496,586 |
0.00, Inf |
>0.9 |
| Low |
372,735,433 |
0.00, Inf |
>0.9 |
# The quantity of interest from a Cox regression model is a hazard ratio (HR). The HR represents the ratio of hazards between two groups at any particular point in time. The HR is interpreted as the instantaneous rate of occurrence of the event of interest in those who are still at risk for the event. It is not a risk, though it is commonly mis-interpreted as such. If you have a regression parameter β, then HR = exp(β).
#
# A HR < 1 indicates reduced hazard of death whereas a HR > 1 indicates an increased hazard of death.
#
# So the HR = 0.59 implies that 0.59 times as many females are dying as males, at any given time. Stated differently, females have a significantly lower hazard of death than males in these data.
### OTHER REPRESENTATION :
s2 <- survfit(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
summary(s2)
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity,
## data = data_survival)
##
## Genetic_Diversity=High
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
##
## Genetic_Diversity=Medium
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 23 2 0.913 0.0588 0.805 1.000
## 5 21 2 0.826 0.0790 0.685 0.996
## 6 19 1 0.783 0.0860 0.631 0.971
##
## Genetic_Diversity=Low
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 34 2 0.941 0.0404 0.865 1.000
## 5 32 4 0.824 0.0654 0.705 0.962
## 6 28 2 0.765 0.0727 0.635 0.921
print(s2)
## Call: survfit(formula = Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity,
## data = data_survival)
##
## n events median 0.95LCL 0.95UCL
## Genetic_Diversity=High 27 0 NA NA NA
## Genetic_Diversity=Medium 23 5 NA NA NA
## Genetic_Diversity=Low 34 8 NA NA NA
# Access to the sort summary table
summary(s2)$table
## records n.max n.start events rmean se(rmean)
## Genetic_Diversity=High 27 27 27 0 6.000000 0.00000000
## Genetic_Diversity=Medium 23 23 23 5 5.739130 0.12627260
## Genetic_Diversity=Low 34 34 34 8 5.764706 0.09355367
## median 0.95LCL 0.95UCL
## Genetic_Diversity=High NA NA NA
## Genetic_Diversity=Medium NA NA NA
## Genetic_Diversity=Low NA NA NA
d <- data.frame(time = s2$time,
n.risk = s2$n.risk,
n.event = s2$n.event,
n.censor = s2$n.censor,
surv = s2$surv,
upper = s2$upper,
lower = s2$lower)
head(d)
## time n.risk n.event n.censor surv upper lower
## 1 6 27 0 27 1.0000000 1.0000000 1.0000000
## 2 4 23 2 0 0.9130435 1.0000000 0.8048548
## 3 5 21 2 0 0.8260870 0.9964666 0.6848395
## 4 6 19 1 18 0.7826087 0.9707088 0.6309579
## 5 4 34 2 0 0.9411765 1.0000000 0.8653187
## 6 5 32 4 0 0.8235294 0.9621763 0.7048611
plot2 <- survminer::ggsurvplot(s2,
pval = TRUE, # show p-value of log-rank test.
#pval.coord = c(1,0.2),
pval.size = 5,
conf.int = TRUE, # show confidence intervals for point estimaes of survival curves.
# conf.int.style = "step", # customize style of confidence intervals
conf.int.style = "ribbon", # customize style of confidence intervals
conf.int.alpha = 0.3, # customize style of confidence intervals
censor = TRUE,
censor.shape = "x",
censor.size = 6,
xlab = "Generation", # customize X axis label.
break.time.by = 1, # break X axis in time intervals by 1
ggtheme = theme_light(), # customize plot and risk table with a theme.
#risk.table = TRUE, # Add risk table
#linetype = "strata", # Change line type by groups
risk.table = "abs_pct", # absolute number and percentage at risk.
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.y.text = FALSE,# show bars instead of names in text annotations in legend of risk table.
risk.table.col = "strata", # Change risk table color by groups
legend.labs = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
legend = c("right"),
legend.title = c("Demographic\nhistory:"),
palette = c("#00AFBB", "#E7B800", "#FC4E07"))
plot2
#
# cowplot::save_plot(file = here::here("figures","Extinction_Survival_analysis_V2.pdf"),
# plot2, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)
# Fit a Cox proportional hazards model
# fit.coxph <- coxph(Surv(Gen_eggs_extinct, status) ~ Genetic_Diversity, data = data_survival)
# ggforest(fit.coxph, data = data_survival)
Population size over
time
Polynomial
analysis
Prelim: test
function
#######################
####################### MODEL INCLUDING GENETIC DIVERSITY WITHOUT GEN 1
#######################
###################### REMOVE GEN 1 ##########################33
#If we dont consider Gen=1
PLOT_Pop_size_average
fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",])
fit <- glm(Nb_adults ~ Generation_Eggs,
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit6 <- glm(Nb_adults~exp(Generation_Eggs),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit7 <- glm(Nb_adults~log(Generation_Eggs),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
#generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(1,6, length=100)
# plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
# log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
# lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
# lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
# lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
# lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
# lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
# lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
data[data$Generation_Eggs>=2,]$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")
# Check the r squared
rsq::rsq(fit,adj=TRUE)
## [1] 0.430959
rsq::rsq(fit2,adj=TRUE)
## [1] 0.5735495
rsq::rsq(fit3,adj=TRUE)
## [1] 0.5723284
rsq::rsq(fit4,adj=TRUE)
## [1] 0.5711286
rsq::rsq(fit5,adj=TRUE)
## [1] 0.5711286
rsq::rsq(fit6,adj=TRUE)
## [1] 0.1285168
rsq::rsq(fit7,adj=TRUE)
## [1] 0.504535
# Best model is fit 2
#data$gen_square <- data$Generation_Eggs^2
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
fit2 <- glm(Nb_adults~Generation_Eggs + gen_square,
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
#Add the interaction with Genetic diversity
fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
gen_square*Genetic_Diversity,
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",],
family = "poisson")
summary(fit2)
##
## Call:
## glm(formula = Nb_adults ~ Generation_Eggs * Genetic_Diversity +
## gen_square * Genetic_Diversity, family = "poisson", data = data[data$Generation_Eggs >=
## 2 & data$Extinction_finalstatus == "no", ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -20.073 -4.181 -0.241 3.315 17.561
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.471305 0.053931 157.076 < 2e-16
## Generation_Eggs -1.688615 0.031961 -52.833 < 2e-16
## Genetic_DiversityMedium -0.593773 0.093819 -6.329 2.47e-10
## Genetic_DiversityLow 0.157769 0.085778 1.839 0.06588
## gen_square 0.184959 0.004076 45.383 < 2e-16
## Generation_Eggs:Genetic_DiversityMedium 0.286248 0.056055 5.107 3.28e-07
## Generation_Eggs:Genetic_DiversityLow -0.241226 0.051392 -4.694 2.68e-06
## Genetic_DiversityMedium:gen_square -0.050817 0.007248 -7.011 2.36e-12
## Genetic_DiversityLow:gen_square 0.019919 0.006595 3.020 0.00253
##
## (Intercept) ***
## Generation_Eggs ***
## Genetic_DiversityMedium ***
## Genetic_DiversityLow .
## gen_square ***
## Generation_Eggs:Genetic_DiversityMedium ***
## Generation_Eggs:Genetic_DiversityLow ***
## Genetic_DiversityMedium:gen_square ***
## Genetic_DiversityLow:gen_square **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 27751 on 352 degrees of freedom
## Residual deviance: 11836 on 344 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: 14121
##
## Number of Fisher Scoring iterations: 5
####### Next step: Final step: add genetic diversity
Analysis
#Test oversdispersion
mod0 <- glm(Nb_adults~Generation_Eggs*Genetic_Diversity +
gen_square*Genetic_Diversity +
Block ,
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",],
family = "poisson")
sqrt(deviance(mod0)/(nobs(mod0)-length(coef(mod0))))
## [1] 5.840179
sigma(mod0)
## [1] 5.840179
# ccl: overdispersion
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs) ,
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",],
family = "poisson")
blmeco::dispersion_glmer(mod1)
## [1] 1.019726
sigma(mod1)
## [1] 1
#Convergence issue
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.0111561
mod1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## Genetic_Diversity + Block + (1 | obs)
## Data: data[data$Generation_Eggs >= 2 & data$Extinction_finalstatus ==
## "no", ]
## AIC BIC logLik deviance df.resid
## 4007.501 4053.899 -1991.751 3983.501 341
## Random effects:
## Groups Name Std.Dev.
## obs (Intercept) 0.6643
## Number of obs: 353, groups: obs, 353
## Fixed Effects:
## (Intercept)
## 9.22506
## Generation_Eggs
## -2.10047
## Genetic_DiversityMedium
## -0.65107
## Genetic_DiversityLow
## 0.10736
## gen_square
## 0.23715
## BlockBlock4
## -0.12145
## BlockBlock5
## -0.48597
## Generation_Eggs:Genetic_DiversityMedium
## 0.32468
## Generation_Eggs:Genetic_DiversityLow
## -0.24262
## Genetic_DiversityMedium:gen_square
## -0.06201
## Genetic_DiversityLow:gen_square
## 0.01594
## convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#Test add Population random
mod0 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs) + (1|Population),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",],
family = "poisson")
max(abs(with(mod0@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.006369406
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.0111561
anova(mod0, mod1, test ="Chisq")
## Data: data[data$Generation_Eggs >= 2 & data$Extinction_finalstatus == ...
## Models:
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## mod1: Genetic_Diversity + Block + (1 | obs)
## mod0: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## mod0: Genetic_Diversity + Block + (1 | obs) + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 12 4007.5 4053.9 -1991.8 3983.5
## mod0 13 3993.9 4044.2 -1984.0 3967.9 15.611 1 7.782e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
########## LRT and Final model:
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs) + (1|Population),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + Genetic_Diversity +
Block + (1|obs) + (1|Population),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod1, mod2, test ="Chisq")
## Data: data[data$Generation_Eggs >= 2 & data$Extinction_finalstatus == ...
## Models:
## mod2: Nb_adults ~ Generation_Eggs + gen_square + Genetic_Diversity +
## mod2: Block + (1 | obs) + (1 | Population)
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## mod1: Genetic_Diversity + Block + (1 | obs) + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod2 9 3997.6 4032.4 -1989.8 3979.6
## mod1 13 3993.9 4044.2 -1984.0 3967.9 11.663 4 0.02004 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(2,6, length = 100),each=3),
Block = "Block4",
gen_square = (rep(seq(2,6, length = 100),each=3))^2,
Estimates = NA)
filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)
#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.3) +
geom_point(data = data,
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
emmeans::emmeans(mod1, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.89 0.0763 Inf 4.71 5.08
## Medium 4.42 0.0930 Inf 4.20 4.65
## Low 4.31 0.0801 Inf 4.12 4.51
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.471 0.120 Inf 3.917 0.0003
## High - Low 0.581 0.110 Inf 5.287 <.0001
## Medium - Low 0.109 0.122 Inf 0.899 0.6408
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
GAMMM Model
Basic model
# https://jacolienvanrij.com/Tutorials/GAMM.html
# install.packages("itsadug")
# packageVersion('mgcv')
# packageVersion('itsadug')
data_dyn <- data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]
str(data_dyn)
## 'data.frame': 355 obs. of 33 variables:
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Old_Label : chr "B3 All red mix " "B3 All red mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "7/7 BE B3 | G5 High 1" "6/2 BE B3 | G4 High 1" "4/28 BE B3 | G3 High 1" "BE B3 3/24 | G2 High1" ...
## $ Genetic_Diversity : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Generation_Parents : int 4 3 2 1 5 5 3 2 1 4 ...
## $ Generation_Eggs : int 5 4 3 2 6 6 4 3 2 5 ...
## $ Date_Census : chr "7/7/21" "6/2/21" "4/28/21" "3/24/21" ...
## $ Date_Start_Eggs : chr "7/7/21" "6/2/21" "4/28/21" "3/24/21" ...
## $ Date_End_Eggs : chr "7/8/21" "6/3/21" "4/29/21" "3/25/21" ...
## $ Image_Box1 : chr "DSC_0779" "DSC_0206" "DSC_0585" "DSC_0920" ...
## $ Image_Box2 : chr "DSC_0780" "DSC_0207" "DSC_0586" "DSC_0921" ...
## $ MC_Box1 : num 26.6 85.4 82.3 266.9 100.1 ...
## $ MC_Box2 : num 5.15 90.72 78.01 149.8 74.38 ...
## $ MC_Total_Adults : num 31.8 176.1 160.3 416.7 174.5 ...
## $ HC_Box1 : int 27 79 83 253 64 47 71 60 101 47 ...
## $ HC_Box2 : int 4 96 82 135 58 80 118 113 214 44 ...
## $ HC_Total_Adults : int 31 175 165 388 122 127 189 173 315 91 ...
## $ Nb_parents : num 175 165 388 100 31 91 173 315 100 189 ...
## $ Growth_rate : num 0.177 1.061 0.425 3.88 3.935 ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 31 175 165 388 122 127 189 173 315 91 ...
## $ Lambda : num 0.177 1.061 0.425 3.88 3.935 ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 337 253 169 85 421 447 279 195 111 363 ...
## $ gen_square : num 25 16 9 4 36 36 16 9 4 25 ...
## $ He_start : num 0.999 0.999 0.999 0.999 0.999 ...
## $ He_lost_gen_t : num 0.984 0.997 0.997 0.999 0.996 ...
## $ He_remain : num 0.971 0.986 0.989 0.992 0.967 ...
## $ He_exp : num 0.968 0.968 0.968 0.968 0.968 ...
## $ He_end : num 0.967 0.967 0.967 0.967 0.967 ...
## $ ID_extinction : chr "extant" "extant" "extant" "extant" ...
# 1-dimensional smooth of Time:
#m1 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs), data=data_dyn)
#Meaning the GAM has failed to construct a sensible smooth term for this data. You can limit the maximum df of the smooth using the k parameter:
m1 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, k=3), data=data_dyn)
#number of ‘knots’. This parameter determines the upper bound of the number of underlying base functions being used to build up the curve. Thus, this parameter constraints the wigglyness of a smooth, or - as a metaphor - the number of bowpoints of a curve.
#Note that the model will base the number of base functions (reflected in the edf of the summary) on the data with the setting for k as upper bound. By default, the value of k for s() is around 9, and for te() and ti() 5 per dimension. Importantly, the value of k should be at most one less than the number of unique data points, otherwise it will fit the density of that predictor.
# interaction with non-isotropicterms:
#m2 <- gam(Nb_adults ~ te(Generation_Eggs, Genetic_Diversity), data=data_dyn)
# decompose the same interaction in main effects + interaction:
# (note: include main effects)
# m3 <- gam(Nb_adults ~ s(Generation_Eggs) + s(Genetic_Diversity) +
# + ti(Generation_Eggs, Genetic_Diversity),
# data=data_dyn)
# s() : for modeling a 1-dimensional smooth, or for modeling isotropic interactions (variables are measured in same units and on same scale)
# te(): for modeling 2- or n-dimensional interaction surfaces of variables that are not isotropic (but see info about d parameter below). Includes ‘main’ effects.
# ti(): for modeling 2- or n-dimensional interaction surfaces that do not include the ‘main effects
###########
########### ADD FACTOR EFFECT:
###########
# How to set up the interaction depends on the type of grouping predictor:
# with factor include intercept difference: Group + s(Time, by=Group).
# with ordered factor include intercept difference and reference smooth: Group + s(Time) + s(Time, by=Group).
# with binary predictor include reference smooth: s(Time) + s(Time, by=IsGroupChildren).
m.factor <- mgcv::gam(Nb_adults ~ Genetic_Diversity + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
summary(m.factor)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Nb_adults ~ Genetic_Diversity + s(Generation_Eggs, k = 3) + s(Generation_Eggs,
## by = Genetic_Diversity, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 172.606 5.715 30.205 < 2e-16 ***
## Genetic_DiversityMedium -46.636 9.015 -5.173 3.92e-07 ***
## Genetic_DiversityLow -55.236 8.160 -6.770 5.56e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Generation_Eggs) 1.714 1.746 129.456 < 2e-16 ***
## s(Generation_Eggs):Genetic_DiversityHigh 1.467 1.667 4.035 0.013248 *
## s(Generation_Eggs):Genetic_DiversityMedium 0.750 0.750 14.965 0.000895 ***
## s(Generation_Eggs):Genetic_DiversityLow 0.750 0.750 8.142 0.013939 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 10/11
## R-sq.(adj) = 0.615 Deviance explained = 62.2%
## GCV = 4473.2 Scale est. = 4375.8 n = 353
#The smooth terms summary shows three smooths that are significantly different from 0.
#We cannot conclude from the summary that the two curves are different from each other.
m.factor <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
summary(m.factor)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 3) +
## s(Generation_Eggs, by = Genetic_Diversity, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 181.701 7.897 23.008 < 2e-16 ***
## Genetic_DiversityMedium -47.619 9.023 -5.277 2.33e-07 ***
## Genetic_DiversityLow -58.396 8.238 -7.088 7.77e-12 ***
## BlockBlock4 -6.052 8.108 -0.746 0.4559
## BlockBlock5 -21.941 9.421 -2.329 0.0204 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Generation_Eggs) 1.714 1.746 130.723 < 2e-16 ***
## s(Generation_Eggs):Genetic_DiversityHigh 1.472 1.670 4.118 0.012347 *
## s(Generation_Eggs):Genetic_DiversityMedium 0.750 0.750 15.121 0.000842 ***
## s(Generation_Eggs):Genetic_DiversityLow 0.750 0.750 8.228 0.013454 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 12/13
## R-sq.(adj) = 0.619 Deviance explained = 62.8%
## GCV = 4452.8 Scale est. = 4330.7 n = 353
###########
########### ADD RANDOM EFFECT:
###########
# random intercepts adjust the height of other modelterms with a constant value: s(Subject, bs="re").
# random slopes adjust the slope ofthe trend of a numeric predictor: s(Subject, Time, bs="re").
# random smooths adjust the trend of a numeric predictor in a nonlinear way: s(Time, Subject, bs="fs", m=1).
m.mixed <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=3) +
s(Population, bs="re") , data=data_dyn)
summary(m.mixed)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 3) +
## s(Generation_Eggs, by = Genetic_Diversity, k = 3) + s(Population,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 181.781 8.979 20.245 < 2e-16 ***
## Genetic_DiversityMedium -47.725 10.259 -4.652 4.79e-06 ***
## Genetic_DiversityLow -58.512 9.362 -6.250 1.29e-09 ***
## BlockBlock4 -6.094 9.217 -0.661 0.5089
## BlockBlock5 -21.800 10.707 -2.036 0.0426 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Generation_Eggs) 1.714 1.746 140.096 < 2e-16 ***
## s(Generation_Eggs):Genetic_DiversityHigh 1.492 1.681 4.504 0.008956 **
## s(Generation_Eggs):Genetic_DiversityMedium 0.750 0.750 16.246 0.000547 ***
## s(Generation_Eggs):Genetic_DiversityLow 0.750 0.750 8.868 0.010341 *
## s(Population) 18.456 66.000 0.389 0.032725 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 83/84
## R-sq.(adj) = 0.646 Deviance explained = 67.3%
## GCV = 4377.5 Scale est. = 4028.3 n = 353
Model
selection
###########
########### AIC
###########
m1 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
m2 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=4), data=data_dyn)
m3 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=5), data=data_dyn)
m4 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, by= Genetic_Diversity, k=3), data=data_dyn)
m5 <- mgcv::gam(Nb_adults ~ Block+
s(Generation_Eggs, by= Genetic_Diversity, k=4), data=data_dyn)
m6 <- mgcv::gam(Nb_adults ~ Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=5), data=data_dyn)
m7 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=5), data=data_dyn)
m8 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=3) +
s(Population, bs="re"), data=data_dyn)
m9 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=4) +
s(Population, bs="re"), data=data_dyn)
m10 <- mgcv::gam(Nb_adults ~ s(Generation_Eggs, by= Genetic_Diversity, k=5)+
s(Population, bs="re"),data=data_dyn)
m11 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, by= Genetic_Diversity, k=3)+
s(Population, bs="re"), data=data_dyn)
m12 <- mgcv::gam(Nb_adults ~ Block +
s(Generation_Eggs, by= Genetic_Diversity, k=4) +
s(Population, bs="re"), data=data_dyn)
m13 <- mgcv::gam(Nb_adults ~ Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
m14 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
m15 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
m16 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=4) +
s(Generation_Eggs, by= Genetic_Diversity, k=4) +
s(Population, bs="re"), data=data_dyn)
m17 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=3) +
s(Population, bs="re"), data=data_dyn)
AIC(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17)
## df AIC
## m1 7.936534 4017.558
## m2 10.081634 4011.546
## m3 11.048573 4011.346
## m4 11.944931 3970.451
## m5 12.082285 4013.676
## m6 11.674632 4012.361
## m7 13.980856 3961.220
## m8 46.262085 3975.939
## m9 49.630906 3965.567
## m10 51.228589 3963.659
## m11 30.598171 3962.136
## m12 50.983869 3966.744
## m13 51.608090 3964.162
## m14 35.378282 3950.158
## m15 33.635622 3947.089
## m16 32.497862 3948.671
## m17 29.161907 3961.031
#https://www.seascapemodels.org/rstats/2021/03/27/common-GAM-problems.html
# Real comparison
mod0 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
mod1 <- mgcv::gam(Nb_adults ~ Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
mod2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn)
mod3 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, by= Genetic_Diversity, k=3) +
s(Population, bs="re"), data=data_dyn)
mod4 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, by= Genetic_Diversity, k=4) +
s(Population, bs="re"), data=data_dyn)
mod5 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
mod7 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=3) +
s(Generation_Eggs, by= Genetic_Diversity, k=3) +
s(Population, bs="re"), data=data_dyn)
mod8 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=4) +
s(Generation_Eggs, by= Genetic_Diversity, k=4) +
s(Population, bs="re"), data=data_dyn)
mod9 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
AIC(mod1, mod2, mod3, mod4, mod5, mod7, mod8, mod9)
## df AIC
## mod1 49.81430 3961.188
## mod2 30.46099 3947.366
## mod3 30.59817 3962.136
## mod4 34.56955 3952.491
## mod5 36.35822 3950.893
## mod7 29.16191 3961.031
## mod8 32.49786 3948.671
## mod9 33.63562 3947.089
mod2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn)
mod9 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by= Genetic_Diversity, k=5) +
s(Population, bs="re"), data=data_dyn)
AIC(mod2, mod8, mod9)
## df AIC
## mod2 30.46099 3947.366
## mod8 32.49786 3948.671
## mod9 33.63562 3947.089
visreg::visreg(mod9, xvar = "Generation_Eggs",
by = "Genetic_Diversity", data = data_dyn,
method = "REML")
visreg::visreg(mod2, xvar = "Generation_Eggs",
by = "Genetic_Diversity", data = data_dyn,
method = "REML")
Testing for
significance
# Three important methods to determine the significance of predictors:
## -1 Model-comparison procedure (using function compareML).
## -2 Test statistics of the smooth terms as presented in the summary.
## -3 Visualization of the smooth terms (e.g., difference plots and difference surfaces) but useful only with significative interaction https://jacolienvanrij.com/Tutorials/GAMM.html
###### 1- Model comparison
mod_1 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by = Genetic_Diversity, k=5) +
s(Population, bs="re"), method="ML", data=data_dyn)
mod_2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn, method='ML')
itsadug::compareML(mod_1, mod_2)
## mod_1: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) +
## s(Generation_Eggs, by = Genetic_Diversity, k = 5) + s(Population,
## bs = "re")
##
## mod_2: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) +
## s(Population, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 mod_2 1975.142 8
## 2 mod_1 1974.369 14 0.772 6.000 0.957
##
## AIC difference: 2.51, model mod_2 has lower AIC.
AIC(mod_1, mod_2)
## df AIC
## mod_1 34.90059 3956.820
## mod_2 30.31497 3954.313
mod_2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn, method='ML')
mod_3 <- mgcv::gam(Nb_adults ~ Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn, method='ML')
itsadug::compareML(mod_3, mod_2)
## mod_3: Nb_adults ~ Block + s(Generation_Eggs, k = 5) + s(Population,
## bs = "re")
##
## mod_2: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) +
## s(Population, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 mod_3 1993.015 6
## 2 mod_2 1975.142 8 17.873 2.000 1.729e-08 ***
##
## AIC difference: 10.73, model mod_2 has lower AIC.
mod_final <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn, method='ML')
###### 2- Summary of the model
# With factors such as Group we cannot use the summary to test differences between groups. However, to investigate differences between groups, we could use difference smooths with ordered factors or binary predictors.
data_dyn$OFGroup_diversity <- as.factor(data_dyn$Genetic_Diversity)
# change factor to ordered factor:
data_dyn$OFGroup_diversity <- as.ordered(data_dyn$OFGroup)
# change contrast to treatment coding (difference curves)
contrasts(data_dyn$OFGroup_diversity) <- 'contr.treatment'
# Inspect contrasts:
contrasts(data_dyn$OFGroup_diversity)
## Medium Low
## High 0 0
## Medium 1 0
## Low 0 1
mod_1 <- mgcv::gam(Nb_adults ~ OFGroup_diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by = OFGroup_diversity, k=5) +
s(Population, bs="re"), method="ML", data=data_dyn)
summary(mod_1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Nb_adults ~ OFGroup_diversity + Block + s(Generation_Eggs, k = 5) +
## s(Generation_Eggs, by = OFGroup_diversity, k = 5) + s(Population,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 181.737 8.773 20.716 < 2e-16 ***
## OFGroup_diversityMedium -47.779 10.023 -4.767 2.83e-06 ***
## OFGroup_diversityLow -58.470 9.148 -6.392 5.70e-10 ***
## BlockBlock4 -6.045 9.005 -0.671 0.5025
## BlockBlock5 -21.621 10.462 -2.067 0.0396 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Generation_Eggs) 3.755 3.963 92.191 <2e-16 ***
## s(Generation_Eggs):OFGroup_diversityMedium 1.217 1.398 0.130 0.6960
## s(Generation_Eggs):OFGroup_diversityLow 1.001 1.001 0.226 0.6353
## s(Population) 17.712 66.000 0.387 0.0208 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.656 Deviance explained = 68.3%
## -ML = 1974.8 Scale est. = 3905.5 n = 353
#We confirm: No interaction significative
# The summary now provides the significance of the difference terms: Medium and High does not show a significantly different trend over Time (F=0.130, edf=1.21, p=0.69), Low and High does not showa significantly different trend over Time (F=0.226, edf=1.001, p=0.02).
mod_2 <- mgcv::gam(Nb_adults ~ OFGroup_diversity + Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), method="ML", data=data_dyn)
summary(mod_2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Nb_adults ~ OFGroup_diversity + Block + s(Generation_Eggs, k = 5) +
## s(Population, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 181.743 8.751 20.768 < 2e-16 ***
## OFGroup_diversityMedium -47.779 9.999 -4.778 2.67e-06 ***
## OFGroup_diversityLow -58.483 9.125 -6.409 5.12e-10 ***
## BlockBlock4 -6.060 8.983 -0.675 0.5004
## BlockBlock5 -21.617 10.436 -2.071 0.0391 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Generation_Eggs) 3.752 3.963 147.751 <2e-16 ***
## s(Population) 17.500 66.000 0.383 0.0205 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.656 Deviance explained = 68.1%
## -ML = 1975.1 Scale est. = 3903.6 n = 353
###### 2- Visualization of the smooth terms
# seful only with significative interaction
# Which is not the case here
# But see: https://jacolienvanrij.com/Tutorials/GAMM.html
Plot
data_dyn <- data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]
mod_final <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block +
s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn)
# Plot block effect
mod_block_P <- tidymv::predict_gam(mod_final,
exclude_terms = list("s(Population)"),
values = list(Population = NULL))
ggplot(data = mod_block_P, aes(Generation_Eggs, fit)) +
facet_wrap(~ Block) +
tidymv::geom_smooth_ci(Genetic_Diversity)
# Plot Diversity effect
mod_A_p <- tidymv::predict_gam(mod_final,
exclude_terms = list("s(Population)", "Block"),
values = list(Population = NULL, Block = NULL))
ggplot(data = mod_A_p, aes(Generation_Eggs, fit)) +
tidymv::geom_smooth_ci(Genetic_Diversity)
## Prediction for main text
#population size at the end
as.data.frame(mod_A_p[mod_A_p$Generation_Eggs==6,])
## Genetic_Diversity Block Generation_Eggs Population fit se.fit
## 1 High Block4 6 High1 139.87617 10.27320
## 2 Medium Block4 6 High1 92.07954 11.62504
## 3 Low Block4 6 High1 81.37454 10.46044
#minimum population size
tempdata <- as.data.frame(mod_A_p)
tapply(mod_A_p$fit,mod_A_p$Genetic_Diversity,min)
## High Medium Low
## 113.09324 65.29661 54.59161
tempdata[tempdata$Genetic_Diversity=="High"&tempdata$fit<=114,]
## Genetic_Diversity Block Generation_Eggs Population fit se.fit
## 91 High Block4 4.448980 High1 113.5282 9.489786
## 94 High Block4 4.530612 High1 113.2052 9.470716
## 97 High Block4 4.612245 High1 113.0932 9.515032
## 100 High Block4 4.693878 High1 113.1954 9.615687
## 103 High Block4 4.775510 High1 113.5150 9.753966
## PLOT
plot_GAM_Ucurve <- ggplot(data = mod_A_p, aes(Generation_Eggs, fit)) +
geom_point(data = data[data$Extinction_finalstatus=="no",],
aes(x = Generation_Eggs,
y = Nb_adults,
group = Genetic_Diversity,
# group = Population,
colour = Genetic_Diversity),
position = position_dodge2(0.3),
size = 0.6, alpha = 0.9) +
tidymv::geom_smooth_ci(Genetic_Diversity,
size = 1.5, linetype = 1) +
ylab("Population size") +
xlab("Generation") +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
# values = c("#685338", "#FEA698", "#7BC7AD")) +
# values = c( "#51392C","#DA5837", "#4C9E97")) +
# values = c( "#00AFBB", "#E7B800", "#FC4E07")) +
# values = c("#9AD3CA", "#929FD1", "#E7C7D7")) +
#values = c("#227C9D", "#FFCB77", "#FE6D73")) +
values = c("#00AFBB","#FC4E07","#E7B800")) +
theme_LO
plot_GAM_Ucurve
#
# cowplot::save_plot(file = here::here("figures","PopulationSize_GAM.pdf"),
# plot_GAM_Ucurve, base_height = 10/cm(1),
# base_width = 15/cm(1), dpi = 610)
#
Variance over
time
Residuals
residuals
#data$gen_square <- data$Generation_Eggs^2
data_G2_all <- data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",]
data_G2_all <- data_G2_all[complete.cases(data_G2_all$Nb_adults), ]
Initial_model<- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs),
data = data_G2_all,
family = "poisson")
data_G2_all$resid <- resid(Initial_model)
PLOT_Pop_size_predict
#PLOT RESIDUALS
ggplot2::ggplot(data_G2_all, aes(x = factor(Generation_Eggs),
y = resid,
colour = Genetic_Diversity)) +
geom_boxplot(outlier.color = NA) +
geom_point(position = position_dodge(0.75), size = 0.7, alpha = 0.5) +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#E7B800","#FC4E07")) +
ylab("Resid(pop size)") +
xlab("Generation") +
theme_LO
###########
########### OPTION A: Test random effects
###########
# Test variances with residuals and random effects:
model_withrandomeffect <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs) + (1|Genetic_Diversity:Generation_Eggs) ,
data = data_G2_all,
family = "poisson")
model_withoutrandomeffect <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs) ,
data = data_G2_all,
family = "poisson")
anova(model_withrandomeffect,model_withoutrandomeffect, test = "Chisq")
## Data: data_G2_all
## Models:
## model_withoutrandomeffect: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## model_withoutrandomeffect: Genetic_Diversity + Block + (1 | obs)
## model_withrandomeffect: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## model_withrandomeffect: Genetic_Diversity + Block + (1 | obs) + (1 | Genetic_Diversity:Generation_Eggs)
## npar AIC BIC logLik deviance Chisq Df
## model_withoutrandomeffect 12 4007.5 4053.9 -1991.8 3983.5
## model_withrandomeffect 13 4009.5 4059.8 -1991.8 3983.5 3e-04 1
## Pr(>Chisq)
## model_withoutrandomeffect
## model_withrandomeffect 0.9859
AIC(model_withrandomeffect)
## [1] 4009.501
AIC(model_withoutrandomeffect)
## [1] 4007.501
#Best model: without random effects
summary(model_withoutrandomeffect)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## Genetic_Diversity + Block + (1 | obs)
## Data: data_G2_all
##
## AIC BIC logLik deviance df.resid
## 4007.5 4053.9 -1991.8 3983.5 341
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.42106 -0.07356 0.01527 0.07992 0.27815
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.4413 0.6643
## Number of obs: 353, groups: obs, 353
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.22506 0.51799 17.809 < 2e-16
## Generation_Eggs -2.10047 0.28072 -7.483 7.29e-14
## Genetic_DiversityMedium -0.65107 0.81458 -0.799 0.424
## Genetic_DiversityLow 0.10736 0.73631 0.146 0.884
## gen_square 0.23715 0.03473 6.828 8.59e-12
## BlockBlock4 -0.12145 0.08305 -1.462 0.144
## BlockBlock5 -0.48597 0.09723 -4.998 5.79e-07
## Generation_Eggs:Genetic_DiversityMedium 0.32468 0.44435 0.731 0.465
## Generation_Eggs:Genetic_DiversityLow -0.24262 0.40166 -0.604 0.546
## Genetic_DiversityMedium:gen_square -0.06201 0.05501 -1.127 0.260
## Genetic_DiversityLow:gen_square 0.01594 0.04970 0.321 0.748
##
## (Intercept) ***
## Generation_Eggs ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## gen_square ***
## BlockBlock4
## BlockBlock5 ***
## Generation_Eggs:Genetic_DiversityMedium
## Generation_Eggs:Genetic_DiversityLow
## Genetic_DiversityMedium:gen_square
## Genetic_DiversityLow:gen_square
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnrt_E Gnt_DM Gnt_DL gn_sqr BlckB4 BlckB5 G_E:G_DM G_E:G_DL
## Genrtn_Eggs -0.972
## Gntc_DvrstM -0.629 0.618
## Gntc_DvrstL -0.696 0.683 0.442
## gen_square 0.936 -0.989 -0.595 -0.658
## BlockBlock4 -0.097 0.000 0.011 0.007 0.000
## BlockBlock5 -0.099 0.012 0.013 0.023 -0.012 0.466
## Gnrt_E:G_DM 0.614 -0.632 -0.978 -0.432 0.625 0.001 -0.007
## Gnrt_E:G_DL 0.679 -0.699 -0.432 -0.978 0.691 0.004 -0.003 0.441
## Gntc_DvrM:_ -0.591 0.625 0.942 0.416 -0.631 -0.001 0.007 -0.989 -0.436
## Gntc_DvrL:_ -0.654 0.691 0.416 0.942 -0.699 -0.004 0.004 -0.437 -0.989
## G_DM:_
## Genrtn_Eggs
## Gntc_DvrstM
## Gntc_DvrstL
## gen_square
## BlockBlock4
## BlockBlock5
## Gnrt_E:G_DM
## Gnrt_E:G_DL
## Gntc_DvrM:_
## Gntc_DvrL:_ 0.441
## convergence code: 0
## Model failed to converge with max|grad| = 0.0430026 (tol = 0.002, component 1)
#Note: should be the same as:
model_lm_withrandomeffects <- lme4::lmer(resid ~ (1|Genetic_Diversity:Generation_Eggs), data = data_G2_all)
lmerTest::rand(model_lm_withrandomeffects)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## resid ~ (1 | Genetic_Diversity:Generation_Eggs)
## npar logLik AIC LRT Df
## <none> 3 -67.275 140.55
## (1 | Genetic_Diversity:Generation_Eggs) 2 -67.275 138.55 -1.4211e-13 1
## Pr(>Chisq)
## <none>
## (1 | Genetic_Diversity:Generation_Eggs) 1
###########
########### OPTION B: Test with Levene's test on residuals
###########
# Test variances with Levene's test
tapply(data_G2_all$resid, list(data_G2_all$Generation_Eggs, data_G2_all$Genetic_Diversity),var)
## High Medium Low
## 2 0.017018250 0.009672648 0.002808944
## 3 0.056168990 0.027477499 0.090848301
## 4 0.096321758 0.171608049 0.178076659
## 5 0.042211123 0.164030698 0.183835255
## 6 0.007039084 0.171966465 0.117733007
# Levene’s test on the varianc
# Using leveneTest()
car::leveneTest(resid ~ interaction(Generation_Eggs, Genetic_Diversity), data = data_G2_all)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 14 2.4148 0.003093 **
## 338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# PairedData::levene.var.test(data_G2_all$dat[df$Genetic_Diversity=="A"], data_G2_all$dat[df$Genetic_Diversity=="B"])
#
# library("PairedData")
#
# leveneTest(data_G2_all$resid[data_G2_all$Genetic_Diversity=="High"&
# data_G2_all$Generation_Eggs==2],
# data_G2_all$resid[data_G2_all$Genetic_Diversity=="Low"&
# data_G2_all$Generation_Eggs==2])
#
#
# ?levene.var.test()
#
#
# ?median()
# ?leveneTest()
#
###########
########### OPTION C: Check heteroscedasticity of residuals
###########
#Breusch-Pagan test
#library(lmtest)
#lmtest::bptest(Initial_model)
#p-value is not significant (i.e., >0.05), we do not reject the null hypothesis. Therefore, we assume that the residuals are homoscedastic
#library(skedastic)
#install.packages("skedastic")
#skedastic::white_lm(Initial_model)
###########
########### OPTION D: Test with Levene's test on data
###########
# Test variances with Levene's test
tapply(data_G2_all$resid, list(data_G2_all$Generation_Eggs, data_G2_all$Genetic_Diversity),var)
## High Medium Low
## 2 0.017018250 0.009672648 0.002808944
## 3 0.056168990 0.027477499 0.090848301
## 4 0.096321758 0.171608049 0.178076659
## 5 0.042211123 0.164030698 0.183835255
## 6 0.007039084 0.171966465 0.117733007
# Levene’s test on the varianc
# Using leveneTest()
car::leveneTest(resid ~ interaction(Generation_Eggs, Genetic_Diversity), data = data_G2_all)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 14 2.4148 0.003093 **
## 338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Raw data
var_raw_data <- ggplot2::ggplot(data[data$Extinction_finalstatus=="no",],
aes(x = factor(Generation_Eggs),
y = Nb_adults,
colour = Genetic_Diversity)) +
geom_boxplot(outlier.color = NA) +
geom_point(position = position_dodge(0.75), size = 0.7, alpha = 0.5) +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse","Intermediate bottleneck","Strong bottleneck"),
values = c("#00AFBB","#FC4E07", "#E7B800")) +
ylab("Population size") +
xlab("Generation") +
geom_text(x=4, y=450, label="**", size = 5, color = "Black") +
theme_LO
var_raw_data
#
#
# cowplot::save_plot(file = here::here("figures","Variance_rawdata.pdf"),
# var_raw_data, base_height = 10/cm(1),
# base_width = 15/cm(1), dpi = 610)
###########
########### OPTION D: Test with Levene's test on data
###########
# Test variances with Levene's test
tapply(data[data$Extinction_finalstatus=="no",]$Nb_adults,
list(data[data$Extinction_finalstatus=="no",]$Generation_Eggs,
data[data$Extinction_finalstatus=="no",]$Genetic_Diversity),var, na.rm=TRUE)
## High Medium Low
## 1 0.000 0.000 0.000
## 2 5442.447 8472.330 5546.986
## 3 6555.977 7844.408 6712.525
## 4 6432.400 4085.585 1243.594
## 5 2658.934 2375.036 1858.627
## 6 1815.000 1940.840 1788.358
# Levene’s test on the varianc
# Using leveneTest()
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no",])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 17 9.4055 < 2.2e-16 ***
## 406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==5,])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.4647 0.6303
## 67
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==4,])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 4.9423 0.009952 **
## 67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==6,])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.3492 0.7065
## 68
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==3,])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5668 0.57
## 68
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&data$Generation_Eggs==2,])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 2.0367 0.1383
## 68
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&
data$Generation_Eggs==4&data$Genetic_Diversity!="Medium",])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 10.689 0.001954 **
## 50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&
data$Generation_Eggs==4&data$Genetic_Diversity!="Low",])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.8783 0.1778
## 42
car::leveneTest(Nb_adults ~ interaction(Generation_Eggs, Genetic_Diversity),
data = data[data$Extinction_finalstatus=="no"&
data$Generation_Eggs==4&data$Genetic_Diversity!="High",])
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.1686 0.1483
## 42
tapply(data[data$Extinction_finalstatus=="no"&
data$Generation_Eggs==4,]$Nb_adults,
data[data$Extinction_finalstatus=="no"&
data$Generation_Eggs==4,]$Genetic_Diversity,
var,
na.rm=TRUE)
## High Medium Low
## 6432.400 4085.585 1243.594
Growth rate G6
Distribution
analysis
head(data_phenotyping)
## Population Week Block ID_Rep Divsersity Population.1 Box Initial.N Start
## 1 High1 5-week Block3 52 High 1 1 30 7/8/21
## 2 High1 5-week Block3 53 High 1 2 30 7/8/21
## 3 High13 5-week Block4 129 High 13 1 30 7/15/21
## 4 High13 5-week Block4 130 High 13 2 30 7/15/21
## 5 High13 5-week Block4 131 High 13 3 30 7/15/21
## 6 High14 5-week Block4 132 High 14 1 30 7/15/21
## End Larvae Pupae Adults Lambda Genetic_Diversity Nb_ttx p_adults
## 1 8/12/21 3 9 97 3.2333333 High 109 0.8899083
## 2 8/12/21 1 1 8 0.2666667 High 10 0.8000000
## 3 8/19/21 2 0 84 2.8000000 High 86 0.9767442
## 4 8/19/21 1 2 104 3.4666667 High 107 0.9719626
## 5 8/19/21 2 1 84 2.8000000 High 87 0.9655172
## 6 8/19/21 1 0 83 2.7666667 High 84 0.9880952
## p_pupae p_larvae He_remain He_start He_exp He_end
## 1 0.08256881 0.027522936 0.9986008 0.9986008 0.9679591 0.9666047
## 2 0.10000000 0.100000000 0.9986008 0.9986008 0.9679591 0.9666047
## 3 0.00000000 0.023255814 0.9986008 0.9986008 0.9786327 0.9772634
## 4 0.01869159 0.009345794 0.9986008 0.9986008 0.9786327 0.9772634
## 5 0.01149425 0.022988506 0.9986008 0.9986008 0.9786327 0.9772634
## 6 0.00000000 0.011904762 0.9986008 0.9986008 0.9766045 0.9752381
#Without considering extinct populations
#tapply(data_phenotyping$Lambda,data_phenotyping$Population,length)
hist(data_phenotyping$Lambda, breaks=50)
# #Which distribution
fitdistrplus::descdist(data_phenotyping$Lambda, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 6.2
## median: 2.1
## mean: 2.213082
## estimated sd: 1.058233
## estimated skewness: 0.5497908
## estimated kurtosis: 3.726017
#############################################################
################## Determine distribution ###################
#############################################################
#Test normal
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda, "norm")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
## 1-mle-norm
## "not rejected"
plot(f1)
f1$aic
## [1] 551.898
#Test log-normal
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda+1, "lnorm")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
## 1-mle-lnorm
## "not rejected"
plot(f1)
f1$aic
## [1] 552.5385
#Test exp
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda, "exp")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
## 1-mle-exp
## "rejected"
plot(f1)
f1$aic
## [1] 669.5117
#Test log-normal
f1 <- fitdistrplus::fitdist(data_phenotyping$Lambda+1, "gamma")
g1 <- fitdistrplus::gofstat(f1)
g1$kstest
## 1-mle-gamma
## "not rejected"
plot(f1)
f1$aic
## [1] 544.6686
## Poisson lognormal model
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
#summary(mlog)
# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
#summary(msqrt)
# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
# # Log Normal
# mlogN <- glm(log(Lambda) ~ Genetic_Diversity + Block,
# family = gaussian, data=data_phenotyping)
## Compare AIC
AIC(mlog, msqrt,mnorm)
## df AIC
## mlog 7 120.8153
## msqrt 7 156.0901
## mnorm 7 528.8983
#Square root a better fits to the data
#But we compare different values: Growth rate and Growth rate+1
# ## Simulate data
x.teo.log <- unlist(simulate(mlog))
x.teo.sqrt <- unlist(simulate(msqrt))
x.teo.norm <- unlist(simulate(mnorm))
# ## QQplot to compared Negative binomial and Poisson log normal distributions
qqplot(data_phenotyping$Lambda, x.teo.log,
main="QQ-plot", xlab="Observed data", ylab="Simulated data",
las=1, xlim = c(0, 10), ylim = c(0, 10), col='blue') ## QQ-plot
points(sort(data_phenotyping$Lambda), sort(x.teo.sqrt), pch=1, col='red')
points(sort(data_phenotyping$Lambda), sort(x.teo.norm), pch=1, col='green')
abline(0,1)
# ## Add legend
legend(0, 7, legend=c("Log", "Sqrt", "Normal"),
col=c("blue", "red", "green"),
pch=1, bty="n")
# ## Normal distribution provides a good fit to the data
# QQ plot
# qqnorm(resid(mlog))
# qqline(resid(mlog))
# qqnorm(resid(msqrt))
# qqline(resid(msqrt))
# qqnorm(resid(mnorm))
# qqline(resid(mnorm)) #best fit
#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)
hist(residuals(msqrt),col="yellow",freq=F) #Seems the better fit
hist(residuals(mnorm),col="yellow",freq=F) #Seems the better fit
#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
## 1-mle-norm
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest
## 1-mle-norm
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
## 1-mle-norm
## "not rejected"
## Compare AIC
AIC(mlog, msqrt, mnorm)
## df AIC
## mlog 7 120.8153
## msqrt 7 156.0901
## mnorm 7 528.8983
#SLog a better fits to the data
Aanalysis
#############################################################
################## Analysis ###################
#############################################################
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
summary(mod0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## Data: data_phenotyping
##
## REML criterion at convergence: 514.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3747 -0.6242 -0.0778 0.5409 3.4451
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2321 0.4817
## Residual 0.7480 0.8649
## Number of obs: 186, groups: Population, 57
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.569709 0.193717 13.265
## Genetic_DiversityMedium -0.873366 0.238587 -3.661
## Genetic_DiversityLow -0.700898 0.222333 -3.152
## BlockBlock4 0.220484 0.214487 1.028
## BlockBlock5 -0.003695 0.249586 -0.015
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -0.460
## Gntc_DvrstL -0.587 0.363
## BlockBlock4 -0.636 0.079 0.175
## BlockBlock5 -0.592 0.089 0.232 0.451
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 5 532.53 548.66 -261.26 522.53
## mod0 7 520.90 543.48 -253.45 506.90 15.635 2 0.0004027 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 2.64 0.135 42.2 2.31 2.98
## Medium 1.77 0.198 51.5 1.28 2.26
## Low 1.94 0.177 58.2 1.51 2.38
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.873 0.239 48.0 3.653 0.0018
## High - Low 0.701 0.223 52.5 3.143 0.0076
## Medium - Low -0.172 0.261 53.0 -0.660 0.7873
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
############
###### Predict
############
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_phenotyping$Genetic_Diversity),
Block = "Block4")
data_predict_extinct$predict <- predict(mod0,
type = "response",
re.form = NA,
newdata = data_predict_extinct)
SUM_GrowthRate <- Rmisc::summarySE(data_phenotyping,
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
Remaining He
difference
Rmisc::summarySE(data[data$Generation_Eggs==2,],
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)
## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2 Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3 Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320
mfull <- lm(He_remain ~ Genetic_Diversity, data=data[data$Generation_Eggs==2,])
mless <- lm(He_remain ~ 1, data=data[data$Generation_Eggs==2,])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
##
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 81 0.00631
## 2 83 3.14572 -2 -3.1394 20162 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.9918 0.001698 81 0.9877 0.9960
## Medium 0.6743 0.001840 81 0.6698 0.6788
## Low 0.5404 0.001513 81 0.5367 0.5441
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.318 0.00250 81 126.829 <.0001
## High - Low 0.451 0.00227 81 198.470 <.0001
## Medium - Low 0.134 0.00238 81 56.200 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Rmisc::summarySE(data[data$Generation_Eggs==2,],
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)
## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2 Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3 Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320
#Calcul for end:
Rmisc::summarySE(data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",],
measurevar = c("He_end"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)
## Genetic_Diversity N He_end sd se ci
## 1 High 27 0.9663114 0.01931641 0.003717444 0.007641317
## 2 Medium 18 0.6378048 0.03493729 0.008234799 0.017373906
## 3 Low 26 0.5090656 0.03443416 0.006753094 0.013908257
mfull <- lm(He_end ~ Genetic_Diversity, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
mless <- lm(He_end ~ 1, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
##
## Model 1: He_end ~ Genetic_Diversity
## Model 2: He_end ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 68 0.06009
## 2 70 2.97522 -2 -2.9151 1649.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.966 0.00572 68 0.952 0.980
## Medium 0.638 0.00701 68 0.621 0.655
## Low 0.509 0.00583 68 0.495 0.523
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.329 0.00905 68 36.316 <.0001
## High - Low 0.457 0.00817 68 55.978 <.0001
## Medium - Low 0.129 0.00912 68 14.124 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates
Growth rate and
He
#Best model with genetic diversity
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod1 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod2 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod5 <- lme4::lmer(Lambda ~ exp(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
MuMIn::model.sel(mod0, mod1,mod2,mod3,mod4,mod5)
## Model selection table
## (Int) Blc Gnt_Dvr He_end log(He_end) exp(He_end) family df
## mod1 0.8154 + 1.766 gaussian(identity) 6
## mod5 0.3445 + 0.8311 gaussian(identity) 6
## mod3 2.5470 + 1.259 gaussian(identity) 6
## mod0 2.5700 + + gaussian(identity) 7
## mod2 2.0770 + gaussian(identity) 5
## mod4 2.0770 + gaussian(identity) 5
## logLik AICc delta weight
## mod1 -257.506 527.5 0.00 0.408
## mod5 -258.097 528.7 1.18 0.226
## mod3 -258.156 528.8 1.30 0.213
## mod0 -257.449 529.5 2.05 0.147
## mod2 -263.551 537.4 9.95 0.003
## mod4 -263.551 537.4 9.95 0.003
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Population'
#Makes sense: best model log He
x <- seq(1:100)
y <- seq(1001:1100)
plot(log(x),y)
#############################################################
################## Analysis ###################
#############################################################
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_end) + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod4 5 532.53 548.66 -261.26 522.53
## mod3 6 522.65 542.01 -255.33 510.65 11.878 1 0.000568 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ log(He_end) + Block + (1 | Population)
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
##
## REML criterion at convergence: 516.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3851 -0.6642 -0.0835 0.5336 3.4375
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2567 0.5067
## Residual 0.7473 0.8644
## Number of obs: 186, groups: Population, 57
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.547160 0.201530 12.639
## log(He_end) 1.259386 0.351126 3.587
## BlockBlock4 0.229875 0.219317 1.048
## BlockBlock5 0.004275 0.254183 0.017
##
## Correlation of Fixed Effects:
## (Intr) lg(H_) BlckB4
## log(He_end) 0.655
## BlockBlock4 -0.625 -0.154
## BlockBlock5 -0.581 -0.197 0.446
#############################################################
################## Predict ###################
#############################################################
data_predict_lambda <- data.frame(He_end =
seq(min(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end),
max(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end),
0.01),
Block = "Block4")
data_predict_lambda$lambda_predict <- predict(mod3,
type = "response",
re.form = NA,
newdata = data_predict_lambda)
PLOT_He_expect <- PLOT_He_Final + geom_line(data = data_predict_lambda,
aes(x = He_end, y = lambda_predict),
linetype = "longdash", col = "grey40", size = 1.25)
#
# #
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"), PLOT_He_expect,
# base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)
#
Adults G6
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 1818.1 1834.3 -904.1 1808.1 181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23572 -0.14532 0.02672 0.14825 0.41394
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.1739 0.4170
## Population (Intercept) 0.1085 0.3294
## Number of obs: 186, groups: ID_Rep:Population, 186; Population, 57
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2943 0.0812 52.884 < 2e-16 ***
## Genetic_DiversityMedium -0.5132 0.1438 -3.569 0.000358 ***
## Genetic_DiversityLow -0.3047 0.1295 -2.352 0.018660 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.563
## Gntc_DvrstL -0.625 0.356
# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population/ID_Rep),
# family = "poisson", data = data_5week)
#dispersion_lme4::glmer(m0) #dispersion ok
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 1826.3 1836.0 -910.15 1820.3
## m0 5 1818.1 1834.3 -904.07 1808.1 12.164 2 0.002284 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 1818.1 1834.3 -904.1 1808.1 181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23572 -0.14532 0.02672 0.14825 0.41394
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.1739 0.4170
## Population (Intercept) 0.1085 0.3294
## Number of obs: 186, groups: ID_Rep:Population, 186; Population, 57
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2943 0.0812 52.884 < 2e-16 ***
## Genetic_DiversityMedium -0.5132 0.1438 -3.569 0.000358 ***
## Genetic_DiversityLow -0.3047 0.1295 -2.352 0.018660 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.563
## Gntc_DvrstL -0.625 0.356
############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.29 0.0812 Inf 4.10 4.49
## Medium 3.78 0.1188 Inf 3.50 4.06
## Low 3.99 0.1011 Inf 3.75 4.23
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.513 0.144 Inf 3.569 0.0010
## High - Low 0.305 0.130 Inf 2.352 0.0489
## Medium - Low -0.208 0.156 Inf -1.340 0.3728
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq")
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 4 1823.6 1836.5 -907.81 1815.6
## m0 6 1820.1 1839.5 -904.07 1808.1 7.4783 2 0.02377 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq")
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 1844.6 1854.3 -919.30 1838.6
## m0 5 1838.5 1854.6 -914.26 1828.5 10.086 2 0.006455 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Adults G2
m0 <- glm(Nb_adults ~ Genetic_Diversity + Block, family = "poisson", data = data[data$Generation_Eggs==2,])
summary(m0)
##
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson",
## data = data[data$Generation_Eggs == 2, ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -18.7057 -2.7621 0.2126 3.3983 11.4454
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.73890 0.01497 383.312 < 2e-16 ***
## Genetic_DiversityMedium -0.19408 0.01628 -11.920 < 2e-16 ***
## Genetic_DiversityLow -0.19278 0.01460 -13.205 < 2e-16 ***
## BlockBlock4 0.10767 0.01579 6.817 9.3e-12 ***
## BlockBlock5 0.17743 0.01600 11.088 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2379.6 on 83 degrees of freedom
## Residual deviance: 2027.9 on 79 degrees of freedom
## AIC: 2667.2
##
## Number of Fisher Scoring iterations: 4
m1<- glm(Nb_adults ~ Block, family = "poisson", data = data[data$Generation_Eggs==2,])
anova(m0,m1,test="Chisq") # difference
## Analysis of Deviance Table
##
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 79 2027.9
## 2 81 2241.4 -2 -213.52 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 5.834 0.01051 Inf 5.809 5.859
## Medium 5.640 0.01243 Inf 5.610 5.670
## Low 5.641 0.01022 Inf 5.617 5.666
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.1941 0.0163 Inf 11.920 <.0001
## High - Low 0.1928 0.0146 Inf 13.205 <.0001
## Medium - Low -0.0013 0.0161 Inf -0.081 0.9964
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
Proportion life
stage
##### P_adults
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 723.41 733.18 -357.70 715.41
## mod0 6 722.93 737.58 -355.46 710.93 4.4827 2 0.1063
rm(mod0,mod1)
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 893.45 906.31 -442.73 885.45
## mod0 6 896.17 915.46 -442.08 884.17 1.284 2 0.5262
rm(mod0,mod1)
##### P_adults
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 692.89 702.67 -342.45 684.89
## mod0 6 692.61 707.26 -340.30 680.61 4.2861 2 0.1173
rm(mod0,mod1)
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 539.89 552.75 -265.95 531.89
## mod0 6 542.26 561.54 -265.13 530.26 1.639 2 0.4407
rm(mod0,mod1)
##### P_larvae
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference and convergence issue
## Data: data_phenotyping_4week
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 441.93 451.71 -216.97 433.93
## mod0 6 445.08 459.74 -216.54 433.08 0.853 2 0.6528
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference and convergence issue
## Data: data_phenotyping
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 737.53 750.39 -364.77 729.53
## mod0 6 740.62 759.91 -364.31 728.62 0.9082 2 0.635
Multinomial life
stage
# data(alligator)
# head(alligator)
# fit <- multinom(food ~ lake+size+sex, data=alligator, weights=count)
# summary(fit$fitted.values)
data_TEMP <- data_phenotyping_all[,c(1:5,11:13,23)]
data_multinom <- tidyr::gather(data_TEMP,"stage","count",6:8)
data_multinom$Divsersity <- as.factor(data_multinom$Divsersity )
head(data_multinom)
## Population Week Block ID_Rep Divsersity He_end stage count
## 1 High1 5-week Block3 52 High 0.9666047 Larvae 3
## 2 High1 5-week Block3 53 High 0.9666047 Larvae 1
## 3 High13 5-week Block4 129 High 0.9772634 Larvae 2
## 4 High13 5-week Block4 130 High 0.9772634 Larvae 1
## 5 High13 5-week Block4 131 High 0.9772634 Larvae 2
## 6 High14 5-week Block4 132 High 0.9752381 Larvae 1
fit <- nnet::multinom(stage ~ Divsersity + Week +
ID_Rep + Population, data=data_multinom, weights=count)
## # weights: 891 (592 variable)
## initial value 21177.949089
## iter 10 value 7279.570688
## iter 20 value 6892.596694
## iter 30 value 6731.300161
## iter 40 value 6673.056095
## iter 50 value 6650.997296
## iter 60 value 6643.800864
## iter 70 value 6640.322734
## iter 80 value 6639.715370
## iter 90 value 6639.359126
## iter 100 value 6639.242868
## final value 6639.242868
## stopped after 100 iterations
summary(fit$fitted.values)
## Adults Larvae Pupae
## Min. :0.1867 Min. :0.0000002 Min. :0.0000015
## 1st Qu.:0.8151 1st Qu.:0.0000715 1st Qu.:0.0025340
## Median :0.9594 Median :0.0158610 Median :0.0209392
## Mean :0.8676 Mean :0.0320594 Mean :0.1003868
## 3rd Qu.:0.9862 3rd Qu.:0.0397470 3rd Qu.:0.1375707
## Max. :1.0000 Max. :0.5357963 Max. :0.6815192
fit
## Call:
## nnet::multinom(formula = stage ~ Divsersity + Week + ID_Rep +
## Population, data = data_multinom, weights = count)
##
## Coefficients:
## (Intercept) DivsersityLow DivsersityMed Week5-week ID_Rep2 ID_Rep3
## Larvae -3.901052 -0.5799254 -0.7792861 -1.384978 -5.682784 -7.3859508
## Pupae -3.081463 1.3120602 -0.1234234 -3.285186 -5.216557 0.3881476
## ID_Rep4 ID_Rep5 ID_Rep6 ID_Rep7 ID_Rep8 ID_Rep9 ID_Rep10
## Larvae 0.2350253 -6.6190065 -3.365056 2.708014 2.686734 -3.484783 2.577294
## Pupae 0.7369804 -0.8498539 -3.616537 2.857288 -5.174837 3.740474 -4.328158
## ID_Rep11 ID_Rep12 ID_Rep13 ID_Rep14 ID_Rep15 ID_Rep16 ID_Rep17
## Larvae -4.906814 1.935461 -0.3904814 0.488392 0.5753569 0.8743354 0.8409919
## Pupae 2.175923 1.432201 -4.1462042 2.185990 2.8109435 2.2998055 -3.5277235
## ID_Rep18 ID_Rep19 ID_Rep20 ID_Rep21 ID_Rep22 ID_Rep23 ID_Rep24
## Larvae 4.432357 -3.141847 -3.856139 -5.4061672 1.60267062 0.5442378 4.3972999
## Pupae 2.739907 -4.891227 1.152743 0.2279343 -0.05506447 0.2597319 0.8666603
## ID_Rep25 ID_Rep26 ID_Rep27 ID_Rep28 ID_Rep29 ID_Rep30
## Larvae 0.561457 0.6544146 0.5416441 -0.3651511 -2.7870179 -2.354927
## Pupae -3.405475 1.8787020 -3.4360034 2.2023796 0.8572163 1.151144
## ID_Rep31 ID_Rep32 ID_Rep33 ID_Rep34 ID_Rep35 ID_Rep36 ID_Rep37
## Larvae 0.6106844 0.5712835 0 1.351590 1.406744 -3.243125 2.621631
## Pupae -0.6697873 1.6316361 0 1.474952 1.697649 -3.346486 -3.858646
## ID_Rep38 ID_Rep39 ID_Rep40 ID_Rep41 ID_Rep42 ID_Rep43 ID_Rep44
## Larvae -4.014978 2.259835 1.052357 0 -3.639545 4.6663160 -3.5701822
## Pupae -2.833921 4.710593 1.581786 0 1.145102 0.5863376 0.1371215
## ID_Rep45 ID_Rep46 ID_Rep47 ID_Rep48 ID_Rep49 ID_Rep50 ID_Rep51
## Larvae 0.6274724 0.7946104 0.5075467 1.285899 -5.863679 1.473193 3.639152
## Pupae 0.6298941 1.3636803 1.2100284 -4.914068 1.217348 1.579528 2.659804
## ID_Rep52 ID_Rep53 ID_Rep54 ID_Rep55 ID_Rep56 ID_Rep57 ID_Rep58
## Larvae 2.677702 3.207160 1.177968 -0.7944315 0.7145318 3.079617 -4.036752
## Pupae 3.916121 4.288607 -4.890511 3.5929641 -4.1264456 4.227368 -2.512159
## ID_Rep59 ID_Rep60 ID_Rep61 ID_Rep62 ID_Rep63 ID_Rep64 ID_Rep65
## Larvae -2.490028 4.796202 -4.641034 -1.929913 -1.319218 3.733072 -4.788136
## Pupae -4.000409 4.248506 4.284235 -3.669149 -2.693820 3.298505 -4.915121
## ID_Rep66 ID_Rep67 ID_Rep68 ID_Rep69 ID_Rep70 ID_Rep71 ID_Rep72
## Larvae -7.017330 3.060861 2.359004 -3.930840 4.427541 0.4542569 -0.1691199
## Pupae 2.347251 -4.708550 2.426601 -3.853992 4.331383 2.3015935 -4.7909903
## ID_Rep73 ID_Rep74 ID_Rep75 ID_Rep76 ID_Rep77 ID_Rep78 ID_Rep79
## Larvae 1.063646 -5.046106 3.907535 -4.577392 3.981346 -5.248417 4.481877
## Pupae 2.367984 1.948639 2.739222 -5.107275 2.539238 -5.606233 3.281810
## ID_Rep80 ID_Rep81 ID_Rep82 ID_Rep83 ID_Rep84 ID_Rep85 ID_Rep86
## Larvae 0 0 0 0 0 1.4398417 0.6681771
## Pupae 0 0 0 0 0 0.6688119 -0.2994476
## ID_Rep87 ID_Rep88 ID_Rep89 ID_Rep90 ID_Rep91 ID_Rep92 ID_Rep93
## Larvae -0.7281906 1.3116933 1.5994132 -0.01299175 0 -4.030158 3.073755
## Pupae 0.8077768 0.2164853 0.4847153 0.69450403 0 1.780968 -0.738383
## ID_Rep94 ID_Rep95 ID_Rep96 ID_Rep97 ID_Rep98 ID_Rep99 ID_Rep100
## Larvae 3.0807121 -1.853249 -3.131841 -0.921824 0 -3.1893152 0
## Pupae 0.4400123 1.702529 1.197579 -3.450526 0 -0.8124809 0
## ID_Rep101 ID_Rep102 ID_Rep103 ID_Rep104 ID_Rep105 ID_Rep106 ID_Rep107
## Larvae 0.7119443 0 1.1894785 -0.3396681 1.5935250 -3.5319853 3.3179488
## Pupae 0.0170997 0 0.4119119 -0.3737761 0.4041902 0.9822233 -0.1501803
## ID_Rep108 ID_Rep109 ID_Rep110 ID_Rep111 ID_Rep112 ID_Rep113 ID_Rep114
## Larvae -3.352196 2.122348 1.429229 2.613800 -3.658613 3.128553 -3.618763
## Pupae 1.119392 0.258575 1.277072 3.006071 -3.696051 2.844143 -3.650847
## ID_Rep115 ID_Rep116 ID_Rep117 ID_Rep118 ID_Rep119 ID_Rep120 ID_Rep121
## Larvae 0 3.2198902 1.8309212 -5.5367772 3.125144 -3.994640 2.713107
## Pupae 0 0.5796406 0.5703267 0.3352911 4.017139 -3.201366 1.211085
## ID_Rep122 ID_Rep123 ID_Rep124 ID_Rep125 ID_Rep126 ID_Rep127 ID_Rep128
## Larvae 2.395045 -4.558337 0 0 0 0 0
## Pupae 3.286414 -6.104341 0 0 0 0 0
## ID_Rep129 ID_Rep130 ID_Rep131 ID_Rep132 ID_Rep133 ID_Rep134 ID_Rep135
## Larvae 0.7688145 -0.1453313 0.1714988 1.433548 2.382642 -4.985808 0.6035421
## Pupae -5.3482942 2.7402744 2.2917514 -4.436682 3.243611 -3.942287 3.6377425
## ID_Rep136 ID_Rep137 ID_Rep138 ID_Rep139 ID_Rep140 ID_Rep141 ID_Rep142
## Larvae 0.3725647 1.470884 0.6168864 0.5868316 -5.408607 2.396927 1.785541
## Pupae -2.1214897 -1.892151 2.4366592 1.0683622 -4.052407 -3.935329 3.721856
## ID_Rep143 ID_Rep144 ID_Rep145 ID_Rep146 ID_Rep147 ID_Rep148
## Larvae 0 2.155239 2.937346 -5.862254 1.264773 -0.3497110
## Pupae 0 2.792150 2.896125 -5.380492 1.775747 -0.5523572
## ID_Rep149 ID_Rep150 ID_Rep151 ID_Rep152 ID_Rep153 ID_Rep154
## Larvae 3.994773e-05 2.221711 0.9995605 -6.340322 1.8621972 1.107054
## Pupae 9.685703e-01 2.204106 1.3875723 -5.911827 0.7904925 2.081951
## ID_Rep155 ID_Rep156 ID_Rep157 ID_Rep158 ID_Rep159 ID_Rep160 ID_Rep161
## Larvae -0.08137172 0.4521213 0.2318772 3.514960 -5.101088 2.125084 -4.753310
## Pupae 0.46348765 0.3775006 1.2331418 3.898539 -4.614923 -3.906614 -4.030408
## ID_Rep162 ID_Rep163 ID_Rep164 ID_Rep165 ID_Rep166 ID_Rep167 ID_Rep168
## Larvae 0.1794121 1.494904 0.01652801 1.472529 0.6700742 0 0
## Pupae 2.8047663 3.561758 1.18116254 1.137290 0.8163059 0 0
## ID_Rep169 ID_Rep170 ID_Rep171 ID_Rep172 ID_Rep173 ID_Rep174 ID_Rep175
## Larvae 4.347524 -5.196228 3.467400 -3.936655 1.4687217 -2.277264 0
## Pupae -4.482572 1.054043 1.092862 1.824772 0.3621998 -3.637726 0
## ID_Rep176 ID_Rep177 ID_Rep178 ID_Rep179 ID_Rep180 ID_Rep181 ID_Rep182
## Larvae 0 -0.5329391 0 2.380801 1.343078 2.851400 -3.543874
## Pupae 0 0.9809249 0 1.840768 -2.838426 -2.567192 -1.627063
## ID_Rep183 ID_Rep184 ID_Rep185 ID_Rep186 ID_Rep187 ID_Rep188 ID_Rep189
## Larvae 1.211254 0.9663694 -0.09957345 0.1184977 -0.7119973 0 2.622016
## Pupae 2.280281 2.5652594 0.96574798 0.9722470 -5.9704332 0 1.066848
## ID_Rep190 ID_Rep191 ID_Rep192 ID_Rep193 ID_Rep194 ID_Rep195 ID_Rep196
## Larvae 0.03957873 -0.9873092 0.3078449 1.575695 1.351085 1.892239 2.035139
## Pupae 0.51743890 0.3698873 0.4552889 4.227675 -3.951720 -3.432143 5.789872
## ID_Rep197 ID_Rep198 ID_Rep199 ID_Rep200 ID_Rep201 ID_Rep202 ID_Rep203
## Larvae -4.458248 -3.326158 0 0 0 0 0.6417168
## Pupae -2.983207 -2.146541 0 0 0 0 3.4166458
## ID_Rep204 ID_Rep205 ID_Rep206 ID_Rep207 ID_Rep208 ID_Rep209 ID_Rep210
## Larvae 0.1886659 -0.01351849 1.204099 2.441037 -4.598971 1.917261 1.856600
## Pupae -4.6919319 3.39713616 -3.849049 2.676855 -4.870301 2.147695 2.466873
## ID_Rep211 ID_Rep212 ID_Rep213 ID_Rep214 ID_Rep215 ID_Rep216 ID_Rep217
## Larvae -3.383011 -4.931778 -2.469360 4.812594 4.730118 -2.443544 -3.760996
## Pupae -3.343382 1.223512 -2.708439 2.529611 -3.825554 -2.668193 3.267193
## PopulationHigh13 PopulationHigh14 PopulationHigh15 PopulationHigh16
## Larvae 0.7949820 -0.5660759 1.843449 -0.02242095
## Pupae -0.3162685 -1.4976160 -4.013641 -0.76085885
## PopulationHigh17 PopulationHigh18 PopulationHigh19 PopulationHigh2
## Larvae 0 -0.7696701 0.9151014 0.3835367
## Pupae 0 0.3077832 2.1919598 -1.2975464
## PopulationHigh20 PopulationHigh21 PopulationHigh22 PopulationHigh23
## Larvae -0.1497994 0.6026268 -1.586127 -0.9539094
## Pupae 0.5522937 2.0741300 -0.716384 -1.5704978
## PopulationHigh24 PopulationHigh25 PopulationHigh26 PopulationHigh27
## Larvae 1.489057 1.4845509 0 1.982130
## Pupae 2.318452 0.8131021 0 2.409463
## PopulationHigh28 PopulationHigh29 PopulationHigh3 PopulationHigh32
## Larvae -0.930248 0 -0.2426038 0.8303826
## Pupae -2.496063 0 -2.4112366 -1.2752861
## PopulationHigh33 PopulationHigh34 PopulationHigh35 PopulationHigh4
## Larvae 1.1905806 -1.7670836 -4.062966 -5.583991
## Pupae -0.4519129 -0.9222608 -2.181869 -1.830637
## PopulationHigh5 PopulationHigh6 PopulationHigh7 PopulationHigh9
## Larvae -2.155828 1.3487831 -2.5011576 0
## Pupae -1.073923 -0.1214129 -0.2045978 0
## PopulationLow1 PopulationLow10 PopulationLow11 PopulationLow12
## Larvae 1.593691 -2.7870179 -2.354927 1.1819679
## Pupae 0.272817 0.8572163 1.151144 0.9618488
## PopulationLow13 PopulationLow14 PopulationLow15 PopulationLow16
## Larvae 1.4398417 1.2516798 1.586421 0
## Pupae 0.6688119 0.7248145 1.179219 0
## PopulationLow17 PopulationLow18 PopulationLow19 PopulationLow2
## Larvae -0.9564029 -2.8262022 0 -1.455092
## Pupae 1.0425846 -0.1104062 0 -2.193611
## PopulationLow20 PopulationLow21 PopulationLow22 PopulationLow23
## Larvae -3.1893152 0 0.7119443 0
## Pupae -0.8124809 0 0.0170997 0
## PopulationLow24 PopulationLow25 PopulationLow26 PopulationLow29
## Larvae 1.1894785 1.25385697 0.6700742 0
## Pupae 0.4119119 0.03041403 0.8163059 0
## PopulationLow3 PopulationLow33 PopulationLow34 PopulationLow35
## Larvae -0.3940592 0 -1.3179595 1.4687217
## Pupae -0.7200337 0 -0.5108957 0.3621998
## PopulationLow4 PopulationLow5 PopulationLow8 PopulationLow9
## Larvae 2.3885948 -2.5656286 1.138041 1.392365
## Pupae -0.3771888 -0.9985766 1.299262 -2.760397
## PopulationMed1 PopulationMed10 PopulationMed11 PopulationMed12
## Larvae 0 1.429229 -1.535023 0
## Pupae 0 1.277072 -1.496684 0
## PopulationMed13 PopulationMed15 PopulationMed16 PopulationMed18
## Larvae -0.4859658 -0.3196806 0 -2.810203
## Pupae 1.4852584 -0.7910699 0 -2.656801
## PopulationMed19 PopulationMed2 PopulationMed21 PopulationMed23
## Larvae 0 -0.4847912 2.380801 1.343078
## Pupae 0 -0.1738838 1.840768 -2.838426
## PopulationMed24 PopulationMed3 PopulationMed4 PopulationMed5
## Larvae -0.6924742 0.8664883 1.052357 0
## Pupae -4.1942544 -1.9819738 1.581786 0
## PopulationMed6 PopulationMed7 PopulationMed8 PopulationMed9
## Larvae -2.543411 1.422083 1.042111 -1.443884
## Pupae 1.868561 1.993574 1.752641 2.210010
##
## Residual Deviance: 13278.49
## AIC: 14046.49
fit0 <- nnet::multinom(stage ~ Divsersity +
ID_Rep + Population , data=data_multinom, weights=count)
## # weights: 888 (590 variable)
## initial value 21177.949089
## iter 10 value 9564.900133
## iter 20 value 8211.500014
## iter 30 value 7849.429014
## iter 40 value 7795.454749
## iter 50 value 7767.734740
## iter 60 value 7753.488447
## iter 70 value 7749.016045
## iter 80 value 7747.279115
## iter 90 value 7746.838532
## iter 100 value 7746.677660
## final value 7746.677660
## stopped after 100 iterations
anova(fit,fit0)
## Likelihood ratio tests of Multinomial Models
##
## Response: stage
## Model Resid. df Resid. Dev Test Df
## 1 Divsersity + ID_Rep + Population 1250 15493.36
## 2 Divsersity + Week + ID_Rep + Population 1248 13278.49 1 vs 2 2
## LR stat. Pr(Chi)
## 1
## 2 2214.87 0
fit1 <- nnet::multinom(stage ~ Week +
ID_Rep + Population, data=data_multinom, weights=count)
## # weights: 885 (588 variable)
## initial value 21177.949089
## iter 10 value 7302.075461
## iter 20 value 6835.928238
## iter 30 value 6705.271706
## iter 40 value 6660.654690
## iter 50 value 6647.757267
## iter 60 value 6642.480165
## iter 70 value 6640.040038
## iter 80 value 6639.517308
## iter 90 value 6639.289994
## iter 100 value 6639.194019
## final value 6639.194019
## stopped after 100 iterations
anova(fit,fit1)
## Likelihood ratio tests of Multinomial Models
##
## Response: stage
## Model Resid. df Resid. Dev Test Df
## 1 Divsersity + Week + ID_Rep + Population 1248 13278.49
## 2 Week + ID_Rep + Population 1248 13278.39 1 vs 2 0
## LR stat. Pr(Chi)
## 1
## 2 0.09769643 0
### New function to test random effect MCLOGIT
# data_TEMP <- data_phenotyping_all[,c(1:5,17:19,23)]
# data_multinom <- tidyr::gather(data_TEMP,"stage","prop",6:8)
# data_multinom$Divsersity <- as.factor(data_multinom$Divsersity )
# head(data_multinom)
#Other example
library("AER")
library("mlogit")
data("TravelMode", package="AER")
head(TravelMode)
## individual mode choice wait vcost travel gcost income size
## 1 1 air no 69 59 100 70 35 1
## 2 1 train no 34 31 372 71 35 1
## 3 1 bus no 35 25 417 70 35 1
## 4 1 car yes 0 10 180 30 35 1
## 5 2 air no 64 58 68 68 30 2
## 6 2 train no 44 31 354 84 30 2
#Here: Each individual can choose across 4 modes of transports: air train bus or car
#We can determine if this choice is random
#or determined by different variables as: vcost, travel, etc.
TM <- mlogit::mlogit.data(data = TravelMode, #dataset
choice = "choice", #Variable containing the choice: yes or no
shape = "long", #type of dataset
alt.levels = "mode") #Variable contianing the list of choice
SUMMARY ALL
ANALYSIS
######### Proba of extinction
mod0 <- glm(y ~ Genetic_Diversity + Block,
data = data_proba_extinction, family = binomial)
mod1 <- glm(y ~ Block ,
data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 81 58.903
## 2 79 46.833 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High -20.07 1855.743 Inf -4451.10 4410.961
## Medium -1.77 0.650 Inf -3.32 -0.216
## Low -1.56 0.532 Inf -2.83 -0.290
##
## Results are averaged over the levels of: Block
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium -18.300 1855.743 Inf -0.010 0.9999
## High - Low -18.506 1855.743 Inf -0.010 0.9999
## Medium - Low -0.206 0.751 Inf -0.274 0.9594
##
## Results are averaged over the levels of: Block
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
#########Time of extinction
mod1 <- glm(Generation_Eggs ~ Genetic_Diversity + Block,
data = data_timeextinction, family = "poisson")
mod0 <- glm(Generation_Eggs ~ Block,
data = data_timeextinction, family = "poisson")
anova(mod0, mod1, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ Genetic_Diversity + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 10 0.95459
## 2 9 0.74888 1 0.20571 0.6502
#########Dynamic over time
# Test effect
mod_1 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Generation_Eggs, by = Genetic_Diversity, k=5) +
s(Population, bs="re"), method="ML", data=data_dyn)
mod_2 <- mgcv::gam(Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k=5) +
s(Population, bs="re"), data=data_dyn, method='ML')
itsadug::compareML(mod_1, mod_2)
## mod_1: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) +
## s(Generation_Eggs, by = Genetic_Diversity, k = 5) + s(Population,
## bs = "re")
##
## mod_2: Nb_adults ~ Genetic_Diversity + Block + s(Generation_Eggs, k = 5) +
## s(Population, bs = "re")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Difference Df p.value Sig.
## 1 mod_2 1975.142 8
## 2 mod_1 1974.369 14 0.772 6.000 0.957
##
## AIC difference: 2.51, model mod_2 has lower AIC.
#########Growth rate
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 5 532.53 548.66 -261.26 522.53
## mod0 7 520.90 543.48 -253.45 506.90 15.635 2 0.0004027 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey")
## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 2.64 0.135 42.2 2.31 2.98
## Medium 1.77 0.198 51.5 1.28 2.26
## Low 1.94 0.177 58.2 1.51 2.38
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.873 0.239 48.0 3.653 0.0018
## High - Low 0.701 0.223 52.5 3.143 0.0076
## Medium - Low -0.172 0.261 53.0 -0.660 0.7873
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
#########Growth rate vs. He
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_end) + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod4 5 532.53 548.66 -261.26 522.53
## mod3 6 522.65 542.01 -255.33 510.65 11.878 1 0.000568 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##### P_adults
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 723.41 733.18 -357.70 715.41
## mod0 6 722.93 737.58 -355.46 710.93 4.4827 2 0.1063
mod0 <- lme4::glmer(p_adults ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 893.45 906.31 -442.73 885.45
## mod0 6 896.17 915.46 -442.08 884.17 1.284 2 0.5262
##### P_adults
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 692.89 702.67 -342.45 684.89
## mod0 6 692.61 707.26 -340.30 680.61 4.2861 2 0.1173
mod0 <- lme4::glmer(p_pupae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 539.89 552.75 -265.95 531.89
## mod0 6 542.26 561.54 -265.13 530.26 1.639 2 0.4407
##### P_larvae
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 441.93 451.71 -216.97 433.93
## mod0 6 445.08 459.74 -216.54 433.08 0.853 2 0.6528
mod0 <- lme4::glmer(p_larvae ~ Genetic_Diversity + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 737.53 750.39 -364.77 729.53
## mod0 6 740.62 759.91 -364.31 728.62 0.9082 2 0.635
All analysis with He
start or He
######### Proba of extinction
#Prepare clean dataset
data_temp <- data[data$Generation_Eggs==6,c("Block","Population","He_start","Extinction_finalstatus")]
data_temp$Extinction_finalstatus <- as.factor(data_temp$Extinction_finalstatus)
data_temp$y <- ifelse(data_temp$Extinction_finalstatus == "yes", 1, 0)
mod0 <- glm(y ~ He_start + Block, data = data_temp, family = binomial)
mod1 <- glm(y ~ Block , data = data_temp, family = binomial)
anova(mod1, mod0, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: y ~ Block
## Model 2: y ~ He_start + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 81 58.903
## 2 80 49.259 1 9.6436 0.0019 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod0)
##
## Call:
## glm(formula = y ~ He_start + Block, family = binomial, data = data_temp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3794 -0.3925 -0.3230 -0.1011 2.6932
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4211 2.1167 0.671 0.5020
## He_start -7.3760 3.1714 -2.326 0.0200 *
## BlockBlock4 0.6707 1.2758 0.526 0.5991
## BlockBlock5 3.0484 1.1382 2.678 0.0074 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 72.388 on 83 degrees of freedom
## Residual deviance: 49.259 on 80 degrees of freedom
## AIC: 57.259
##
## Number of Fisher Scoring iterations: 6
#########Time of extinction
mod1 <- glm(Generation_Eggs ~ He_start + Block,
data = data_timeextinction, family = "poisson")
mod0 <- glm(Generation_Eggs ~ Block,
data = data_timeextinction, family = "poisson")
anova(mod0, mod1, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: Generation_Eggs ~ Block
## Model 2: Generation_Eggs ~ He_start + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 10 0.95459
## 2 9 0.76852 1 0.18607 0.6662
summary(mod1)
##
## Call:
## glm(formula = Generation_Eggs ~ He_start + Block, family = "poisson",
## data = data_timeextinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.42430 -0.14557 0.03684 0.07219 0.50431
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4209 1.5248 1.588 0.112
## He_start -0.9243 2.1583 -0.428 0.668
## BlockBlock4 -0.1473 0.5277 -0.279 0.780
## BlockBlock5 -0.3318 0.4810 -0.690 0.490
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1.40752 on 12 degrees of freedom
## Residual deviance: 0.76852 on 9 degrees of freedom
## AIC: 53.687
##
## Number of Fisher Scoring iterations: 4
#########Dynamic over time
# Test effect
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*He_start + gen_square*He_start +
Block + (1|obs) + (1|Population),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + He_start +
Block + (1|obs) + (1|Population),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",], family = "poisson")
#Do not converge
#anova(mod1, mod2, test ="Chisq")
#########Growth rate
mod0 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod1 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod0,mod1,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 5 532.53 548.66 -261.26 522.53
## mod0 6 521.90 541.26 -254.95 509.90 12.627 1 0.0003802 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#########Growth rate vs. He
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_end) + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod4 5 532.53 548.66 -261.26 522.53
## mod3 6 522.65 542.01 -255.33 510.65 11.878 1 0.000568 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##### P_adults
mod0 <- lme4::glmer(p_adults ~ He_end + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 723.41 733.18 -357.70 715.41
## mod0 5 721.94 734.16 -355.97 711.94 3.4655 1 0.06266 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lme4::glmer(p_adults ~ He_end + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_adults ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_adults ~ 1 + Block + (1 | Population)
## mod0: p_adults ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 893.45 906.31 -442.73 885.45
## mod0 5 894.05 910.13 -442.03 884.05 1.4012 1 0.2365
##### P_adults
mod0 <- lme4::glmer(p_pupae ~ He_end + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 692.89 702.67 -342.45 684.89
## mod0 5 691.48 703.70 -340.74 681.48 3.4111 1 0.06476 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lme4::glmer(p_pupae ~ He_end + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_pupae ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_pupae ~ 1 + Block + (1 | Population)
## mod0: p_pupae ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 539.89 552.75 -265.95 531.89
## mod0 5 540.28 556.35 -265.14 530.28 1.6188 1 0.2033
##### P_larvae
mod0 <- lme4::glmer(p_larvae ~ He_end + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population),
data = data_phenotyping_4week,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping_4week
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 441.93 451.71 -216.97 433.93
## mod0 5 443.34 455.55 -216.67 433.34 0.5929 1 0.4413
mod0 <- lme4::glmer(p_larvae ~ He_end + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
mod1 <- lme4::glmer(p_larvae ~ 1 + Block + (1|Population),
data = data_phenotyping,family = binomial(link = "logit"),
weights = Nb_ttx)
anova(mod0,mod1,test="Chisq") # NO difference
## Data: data_phenotyping
## Models:
## mod1: p_larvae ~ 1 + Block + (1 | Population)
## mod0: p_larvae ~ He_end + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 737.53 750.39 -364.77 729.53
## mod0 5 738.90 754.97 -364.45 728.90 0.6313 1 0.4269
